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THE  VIEWS  AND  CONCLUSIONS  CONTAINED  IN  THIS  DOCUMENT  ARE  THOSE  OF  THE 
AUTHORS  AND  SHOULD  NOT  BE  INTERPRETED  AS  REPRESENTING  THE  OFFICIAL  POLICIES, 
EITHER  EXPRESSED  OR  IMPLIED,  OF  THE  U.S.  AIR  FORCE 


JOINT  INVERSION  OF  HYDROLOGIC  AND  GEOPHYSICAL  DATA  FOR 
PERMEABILITY  DISTRIBUTION  IN  AN  ALLUVIAL  AQUIFER 


Warren  Barrash,  Paul  Donaldson,  Paul  Michaels,  Jack  Pelton  and  Scott  Urban 
Center  for  Geophysical  Investigation  of  the  Shallow  Subsurface 
Boise  State  University,  Boise,  Idaho 


Abstract 

We  are  in  the  first  year  of  a  three-year  field  and  modeling  study  to  map  permeability  and 
hydrostratigraphy  in  a  heterogeneous  alluvial  aquifer  by  extending  control  at  wells  with 
non-invasive  geophysical  investigation  methods.  Our  efforts  build  upon  recent  work  [1,2] 
which  used  laboratory -derived  relationships  between  permeability  and  seismic  P-velocity, 
and  joint  inversion  of  synthetic  data  sets,  to  demonstrate  the  feasibility  of  determining  the 
distribution  of  permeability  in  a  heterogeneous  aquifer  without  requiring  information  from 
a  dense  network  of  wells. 

The  aquifer  under  study  is  a  heterogeneous  alluvial  system  with  50-70  ft  (15-20  m)  of  sat¬ 
urated  thickness  above  a  clay  layer  at  the  25-acre  Capital  Station  remediation  site  in  Boise, 
Idaho  (Fig.  1A).  Our  approach  is  to  use  seismic  methods  (reflection,  refraction,  surface 
wave  dispersion),  ground  penetrating  radar,  resistivity,  and  time-domain  electromagnetic 
methods  (Fig.  IB)  to  determine  site-wide  subsurface  distributions  of  seismic  velocity  (P- 
wave  and  S-wave),  dielectric  constant,  and  resistivity.  Profiles  of  hydrostratigraphy  and 
permeability  are  being  obtained  from  wells  (Fig.  1C)  with  geophysical  logging,  sampling/ 
coring,  vertical  seismic  profiling  and  multiple-well  pumping  tests. 

Borehole  and  laboratory  data  will  be  used  to  develop  a  locally  valid  empirical  relationship 
between  permeability  and  geophysical  parameter  space  expressed  in  terms  of  radial  basis 
functions.  An  important  part  of  the  analysis  will  be  determining  which  parts  of  geophysi¬ 
cal  parameter  space  make  the  greatest  contributions  in  the  determination  of  permeability. 
The  empirical  relationship  will  be  tested  by  examining  its  predictive  ability  at  individual 
well  locations,  and  then  it  will  be  applied  to  generate  a  distribution  of  permeability 
throughout  the  study  area. 

Geologic  and  Geophysical  Logs.  Drilling  logs  from  existing  wells  describe  a  generalized 
two-layer  stratigraphy  in  the  saturated  zone  with  cobbles  in  a  sand  matrix  overlying  a  sand 
unit,  and  these  data  suggest  a  directional  depositional  fabric  trending  NNW.  More  detailed 
interpretations  of  stratigraphy  and  physical  properties  at  individual  wells  come  from  bore¬ 
hole  geophysical  logs  (e.g.,  Figs.  2-3),  and  from  samples  or  core  from  new  wells  at  the 
site.  Neutron,  natural  gamma,  resistivity,  conductivity  and  temperature  logs  were  run  in  8 
deep  and  9  shallow  wells  in  1994  by  the  U.S.  Geological  Survey  Borehole  Geophysics 
Research  Project  (Fig.  1C);  a  second  round  of  logging  is  planned  for  1995  to  investigate 
remaining  monitoring  wells  and  new  wells.  Core  samples  (disturbed)  from  sands  in  one 
well  are  being  analyzed  for  Vp,  Vs,  grain  density,  grain  size  distribution,  porosity  and  per¬ 
meability;  we  are  collaborating  with  rock  physics  researchers  at  Stanford  [1,  2]  on  analysis 
and  interpretation  of  these  samples. 
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Aquifer  Tests.  Fourteen  multiple-well  aquifer  tests  have  been  run  at  the  site  where  wells 
are  screened  in  different  intervals  (Fig.  1C),  including  tests  in  four  areas  where  two  or 
more  wells  were  used  alternately  for  pumping  and  observation  (e.g.,  Fig.  2).  Simple  1-  or 
2-layer  analytical  or  numerical  models  of  the  system  do  not  adequately  match  observed 
behavior.  Modeling  of  local  areas  pumped  with  several  wells  using  the  axisymmetric 
option  of  MODFLOW  gives  reasonable  matches  if  additional  layers  are  added  to  delay 
onset  of  drawdown  in  observation  wells.  Where  available,  borehole  geophysical  logs  and 
core  samples  provide  control  to  assist  development  of  layered  models  (Fig.  2B).  Perme¬ 
ability  (or  hydraulic  conductivity)  of  layers  ranges  over  3  orders  of  magnitude  with  clean 
sands  being  the  most  permeable  units.  At  least  7  more  tests  will  be  run  in  remaining  mon¬ 
itoring  wells,  and  additional  tests  will  be  run  as  new  wells  are  emplaced. 

Vertical  Seismic  Profiles.  Vertical  seismic  profiles  of  P-  and  S-wave  velocity  in  three  bore¬ 
holes  (Fig.  1C)  were  run  in  1994  using  two  downhole  geophone  tools.  Repeatable  hori¬ 
zontal  and  vertical  hammer  sources  were  used  to  acquire  3-component  data  from  which 
broad-band  first  arrival  time  measurements  have  been  made  (e.g.,  Fig.  3).  Figure  3C 
shows  a  vertical  travel-time  profile.  As  expected,  the  water  table  provides  the  largest  dis¬ 
continuity  in  P-wave  velocity.  However,  the  similar  change  in  S-wave  velocity  was  unex¬ 
pected  because  S-waves  should  be  relatively  insensitive  to  variations  in  pore  fluid  content. 
In  addition,  the  P-wave  reflection  image  suggests  that  mode  conversions  are  occurring  at 
the  water  table.  As  a  result,  we  are  evaluating  the  potential  value  of  offset-dependent 
observations  for  future  surface  seismic  work.  In  Figure  3B,  the  silty-sand  layer  (15-19  m 
depth)  appears  to  generate  a  good  S-wave  reflection,  making  this  also  a  possible  objective 
for  future  S-wave  reflection  work.  Fine  scale  velocity  anomalies  within  the  saturated  zone 
are  also  apparent  in  the  travel  time  data  (Fig.  3C).  These  fluctuations  are  most  noticeable 
on  the  P-wave  times.  We  hope  to  relate  the  P-wave  variations  to  changes  in  hydrologic 
properties  by  using  the  S-waves  as  a  datum. 

Future  plans  include  development  of  a  downhole  geophone  tool  that  will  fit  in  2-inch  bore¬ 
holes  and  couple  with  a  mechanical,  rather  than  pneumatic,  mechanism  in  order  to  collect 
high-quality  S-wave  data  in  standard  monitoring  wells.  Also  we  plan  to  explore  the  atten¬ 
uation  and  phase-velocity  dispersion  in  the  full  waveform  data.  We  hope  to  achieve  verti¬ 
cal  resolution  of  4  m  from  these  measurements. 

Surface  Wave  Dispersion.  Dispersion  of  seismic  surface  waves  may  be  used  to  infer  near¬ 
surface  physical  properties.  In  particular,  the  profile  of  shear-wave  velocity  (V^)  with 
depth  may  be  inferred  from  the  dispersion  of  Rayleigh  waves.  Measurements  of  Ray¬ 
leigh-wave  dispersion  were  carried  out  at  the  Capital  Station  site  to  independently  confirm 
Vs  data  obtained  from  downhole  seismic  VSP  experiments,  and  to  test  the  feasibility  of 
using  surface-wave  techniques  as  a  low-cost  method  for  extending  VSP  results  away  from 
boreholes. 

The  Capital  Station  seismic  experiments  show  that  Rayleigh  waves  are  easily  generated 
and  recorded  along  profiles  of  approximately  100-m  length  virtually  anywhere  on  the  site. 
Preliminary  interpretations  of  phase-velocity  dispersion  curves  for  the  fundamental  Ray- 
leigh-wave  mode  are  based  on  a  simple  single  layer  over  half-space  velocity  model.  These 
interpretations  show  that  estimates  of  Vs  are  consistent  with  estimates  derived  from  nearby 
VSP  experiments,  but  that  VSP-derived  interface  depths  are  greater  than  Rayleigh-wave 
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interpretations  by  a  factor  of  160-200%  (Table  1).  One  explanation  for  the  differences  in 

TABLE  1.  Comparison  of  shear-wave  velocity  estimates  derived  from  Rayleigh-wave  dispersion 
data  and  VSP  measurements  in  boreholes.  Velocity  model  is  single  layer  over  a  half-space.  Wells 
chosen  for  comparison  are  those  with  VSP  measurements  located  closest  to  seismic  lines  2W  (well 
14)  and  2N  (well  ST3). 


Layer  Vs 

Half-Space  Vs 

Interface  Depth 

Line  2W 

1 80  m/s 

2.8  m 

Well  14 

180  m/s 

4.5  m 

Line  2N 

1 50  m/s 

3.1  m 

Well  ST3 

1 70  m/s 

5.9  m 

interface  depth  may  be  the  fact  that  the  Rayleigh-wave  estimates  represent  averages  over  a 
100-m  profile,  whereas  the  VSP  data  are  point  measurements  offset  to  one  end  of  the  seis¬ 
mic  lines.  However,  because  the  VSP  data  indicate  that  the  interface  is  the  water  table,  and 
because  it  is  unlikely  that  the  average  water  table  depth  beneath  the  seismic  line  differs  sig¬ 
nificantly  from  the  depth  at  the  well,  it  appears  that  a  systematic  bias  may  exist  in  the  pre¬ 
liminary  Rayleigh-wave  interpretations.  It  is  probable  that  this  discrepancy  is  related  to 
non-uniqueness  in  the  inversion  procedure  and  is  under  active  investigation. 

An  example  of  the  Rayleigh-wave  data  and  dispersion  analysis  is  shown  in  Figure  4  where 
seismic  records  for  lines  2W  and  2N  (about  300  m  apart,  Fig.  IB)  are  presented  for  com¬ 
parison.  Both  lines  were  shot  with  the  same  source-receiver  geometry  (Fig.  4A).  Rayleigh 
waves  dominate  the  seismograms  within  the  triangular  windows  bracketed  by  the  white 
lines,  but  the  window  for  line  2N  is  longer  in  duration  than  for  line  2W,  indicating  more 
dispersion  and  a  concomitant  greater  stretching  of  the  waveforms.  This  difference  is 
reflected  in  the  phase-velocity  dispersion  curves  where  the  overall  slope  of  the  2N  data  is 
greater  than  the  slope  of  the  2W  data  (Fig.  4B).  Despite  the  clear  differences  in  the  recorded 
data  for  lines  2W  and  2N,  the  associated  velocity  models  are  not  greatly  different  (Fig.  4C), 
suggesting  that  dispersion  is  a  sensitive  phenomenon  for  the  range  of  model  parameters 
under  consideration.  These  preliminary  results  indicate  that  additional  work  on  the  surface- 
wave  dispersion  method  is  warranted,  with  emphasis  on:  increasing  the  bandwidth  (espe¬ 
cially  at  lower  frequencies);  inclusion  of  higher  modes  beyond  the  fundamental  mode; 
observations  of  horizontally  polarized  surface  waves  (Love  waves);  and  consideration  of 
more  complicated  velocity  models. 

Resistivity  Soundings.  Preliminary  resistivity  data  from  8  profiles  acquired  in  the  fall  of 
1994  (Fig.  IB)  have  been  inverted,  and  compared  and  correlated  with  lithologic  logs  and 
geophysical  logs.  Some  soundings  (in  northwest  of  site)  appear  to  be  contaminated  by  lat¬ 
eral  disturbances.  The  best  correlations  are  observed  in  the  southeast  part  of  the  site  where 
the  soundings  reflect  changes  observed  in  lithologic  and  geophysical  logs.  Sounding  2N 
(Fig.  IB)  has  been  included  as  an  example  (Fig.  5).  These  observations  predict  probable 
success  for  our  planned  Transient  EM  (TEM)  experiments.  The  TEM  technique  is  noted  for 
minimizing  lateral  effects  from  vertical  soundings.  We  will  be  experimenting  on-site,  seek¬ 
ing  the  minimum-diameter  transmitter  loop  feasible.  Minimizing  this  parameter  corre¬ 
sponds  to  optimizing  lateral  resolution  in  our  attempts  to  image  the  subsurface  resistivity 
distribution  which  in  turn  can  be  correlated  to  porosity/permeability. 
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Ground  Penetrating  Radar.  Tests  of  ground  penetrating  radar  (GPR)  along  a  single  profile 
(Fig.  IB)  suggest  that  GPR  measurements  will  produce  little  useful  information.  The  tests 
were  conducted  using  a  Sensors  and  Software  pulseEKKO  IV  unit  with  separable  source 
and  receiver.  Common  midpoint  experiments  using  25-MHz,  50-MHz  and  100-MHz 
antennas  failed  to  provide  convincing  evidence  of  any  subsurface  reflectors.  Independent 
resistivity  measurements  (Fig.  5)  indicate  that  the  unsaturated  zone  is  relatively  conduc¬ 
tive  (5-25xl0'3  S/m),  thereby  limiting  the  usefulness  of  GPR  through  severe  attenuation 
of  the  electromagnetic  waves. 
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Figure  1 .  Capital  Station  site  (within  dashed  line),  downtown  Boise,  Idaho.  A)  Shallow  wells  are 
<25  ft  deep  and  open  across  the  water  table  in  the  cobble  unit.  Deep  wells  are  50-85  ft  deep  and 
open  within  the  sand  unit.  B)  Surface  geophysical  lines  run  in  1994.  C)  Wells  tested 
geophysically  and  hydrologically  in  1994. 
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Symbols:  g  (natural  gamma  log),  n  (neutron  log). 
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Figure  4.  Seismic  surface  wave  data  and  preliminary  interpretations  for  Lines  2W  and  2N  at  Capital  Station  site:  A)  48-channel 
records  with  surface-wave  window  indicated  by  white  lines,  traces  scaled  individually;  B)  phase-velocity  dispersion  curves;  C)  pre¬ 
liminary  single  layer  over  half-space  shear-wave  velocity  models  interpreted  from  phase-velocity  dispersion  curves. 
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Figure  5.  Schlumberger  resitivity  sounding  2N  (see  Fig.  IB)  centered  on  well  SA1.  Depth  scale  for  resistivity  interpretation 
is  logarithmic;  depth  scale  for  lithologic  log  is  linear.  Note  that  lithologic  log  starts  below  a  depth  of  5  m  (below  the  top 
of  the  cobble  layer).  The  resistivity  data  do  not  distinguish  layers  beneath  the  base  of  the  cobbles. 
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Summary  The  adjoint  state  computation  is  essential  for  iterative  linearized  inver¬ 
sion  using  full  wave  equation  simulation.  For  dissipative  systems  like  viscoacoustics 
however  this  computation  appears  to  require  storage  of  the  entire  time  history  of 
the  forward  reference  field,  or  else  its  recomputation  at  every  time  step,  both  pro¬ 
hibitively  expensive  in  a  multidimensional  setting.  The  computational  complexity 
and  storage  requirements  can  both  be  reduced  to  a  level  only  exceeding  that  of  the 
forward  computation  by  a  logarithmic  factor,  through  retention  of  a  few  carefully 
selected  intermediate  time  levels  during  the  forward  calculation.  The  feasibility  of 
the  resulting  adjoint  state  computation  is  illustrated  by  linearized  viscoacoustic  2D 
inversion  of  a  small  part  of  a  field  seismic  data  set. 

The  Rice  Site  Characterization  Project  The  project  plan  includes  two  major 
emphases:  development  of  modeling  methods  for  ground  penetrating  radar  waves, 
and  employment  of  these  methods  to  evaluate  the  potential  of  inversion  of  GPR  data 
as  a  tool  for  hazardous  waste  site  characterization  and  monitoring.  Modeling  issues 
to  be  studied  include  accuracy  and  efficiency  of  explicit  finite  element  methods  in 
smooth  and  irregular  models,  and  approximation  of  complex  boundary  conditions 
and  antenna  simulation.  Inversion  issues  will  echo  those  studied  in  The  Rice  Inver¬ 
sion  Project  for  exploration  scale  seismology:  estimation  of  the  velocity  field  (index 
of  refraction),  calibration  of  sources,  and  effects  of  attenuation  and  dispersion  on 
linearized  inversion  for  wavelength  scale  features.  As  attenuation  and  dispersion  are 
of  considerably  greater  importance  in  shallow  seismology  and  GPR  studies  than  in 
exploration  seismology,  this  presentation  reviews  some  recent  work  or  our  group  on 
viscoelastic  inversion  and  its  application  to  seismic  field  data. 


Linearized  Inversion  The  linear  viscoelastic  wave  equation  may  be  written  in  the 
form 

^  =  LU  +  F-  U,  F  =  0,  t  «  0 
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Here  the  vector  U  consists  of  particle  velocity  and  stress  components,  as  well  as  so 
called  memory  variables  which  permit  a  differential  implementation  of  the  viscous 
stress-strain  relation.  The  resulting  system  is  quite  similar  to  Maxwell’s  equations 
with  Debye-type  relaxation  modeling.  The  coefficients  of  the  partial  differential  op¬ 
erator  L  are  (functions  of)  the  mechanical  parameters  of  the  viscoelastic  medium, 
such  as  density,  wave  velocities,  and  relaxation  times,  and  express  Newton’s  law  and 
the  constitutive  relation  (generalized  Hooke’s  law)  of  the  material.  The  vector  valued 
function  F  models  the  energy  source  mechanism.  Finally  let  M  denote  the  measure¬ 
ment  operator,  i.e.  sampling  at  receivers. 

The  linearized  forward  map  relates  perturbations  in  the  coefficients  of  L  to  the  cor¬ 
responding  perturbations  in  MU .  The  least  squares  formulation  of  the  linearized 
inverse  problem  consists  in  minimizing  the  mean  square  error  \\MSU  —  8D ||2,  where 
8D  is  the  residual  data  not  already  explained  by  the  current  reference  model  (coeffi¬ 
cients  of  L).  Minimization  by  descent  methods  requires  calculation  of  the  gradient, 
which  in  turn  demainds  computation  of  the  adjoint  operator  to  8L  i-+  M8U.  The 
following  trick  computes  the  adjoint:  solve  the  final  value  problem 

dW 

—  =  -L*W  +  M*$ ;  W  =  Q,t  »  0 
at 

Then  integration  by  parts  shows  that  <  $,M8U  >=<  W,8LU  >  where  the  angle 
brackets  represent  appropriate  L 2  inner  products.  From  this  formula  one  reads  off 
the  adjoint  of  the  linearied  map,  hence  the  gradient  of  the  mean  square  error  (using 
$  =  8D). 

For  more  details  on  of  linearized  inversion  in  viscoacoustic  media  see  [1],  also  Tarantola 

[4]- 

Recomputation  scheme  As  one  sees  from  last  two  equations,  calculation  of  the 
adjoint  map  involves  accessing  the  solutions  of  initial  and  final  value  problems  at  the 
same  instant  of  time.  One  could  recompute  the  solution  of  the  initial  value  problem 
for  each  time  step  (after  discretization!),  or  retain  in  core  memory  the  entire  solu¬ 
tion  of  the  initial  value  problem,  not  just  the  simulated  measurements.  For  problems 
in  more  than  one  space  variable,  both  alternatives  are  infeasibly  expensive.  A  third, 
compromise  alternative  saves  a  number  of  “snapshots”  (intermediate  complete  states, 
or  Cauchy  data)  and  restarts  the  computation.  We  call  such  an  approach  a  recom¬ 
putation  scheme.  Griewank  [2]  devised  a  recomputation  scheme  which  is  optimal  in 
both  computation  effort  and  memory  use.  It  allows  for  logarithmic  growth  in  com¬ 
putation  effort  and  memory  with  linearly  increasing  problem  size.  We  have  adapted 
Griewank’s  scheme  to  linearized  viscoacoustic  inversion. 

The  algorithm  is  very  efficient.  For  10,000  timesteps,  memory  available  for  only 
10  snapshots  (out  of  10,  000  for  a  “store  it  all”  algorithm)  is  sufficient  to  limit  the 
computational  effort  to  roughly  3  times  more  than  for  a  “store  it  all”  algorithm.  For 
merely  5  available  snapshots,  roughly  5  times  larger  computational  effort  is  sufficient. 
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Viscoacoustic  Inversion  of  the  Mobil  AVO  data  set  In  1994  Mobil  Research 
and  Development  Corporation  made  a  marine  data  set  publicly  available  for  Am¬ 
plitude  -  Versus  -  Offset  (“AVO”  -  roughly  synonymous  with  “linearized  inversion”) 
studies.  Many  research  groups,  including  TRIP,  presented  initial  results  at  a  Soci¬ 
ety  of  Exploration  Geophysicists  workshop  in  1994  entitled  “Comparison  of  seismic 
inversion  methods  on  a  single  real  dataset”,  organized  by  Robert  Keys  and  Douglas 
Foster.  The  data  is  supplemented  by  measurements  in  several  nearby  wells. 

To  demonstrate  the  viscoacoustic  algorithm  we  performed  two  crude  linearized  inver¬ 
sions  of  a  (small)  part  of  the  data  set,  one  with  physically  realistic  attenuation  and 
one  with  no  attenuation  (elastic  modeling).  We  selected  nine  shot  records  located 
close  to  one  of  the  wells.  Preprocessing  included  multiple  reflection  removal  (to  make 
the  linearization  hypothesis  more  accurate),  amplitude  compensation  (to  enable  2D 
treatment  of  data  recorded  in  3D),  editing,  downfiltering  to  allow  use  of  a  coarser 
finite  difference  grid,  and  source  ghost  removal.  Background  (smooth)  velocity,  bulk 
modulus,  density,  and  attenuation  models  were  derived  from  the  well  log  and  vertical 
seismic  profile  data.  The  background  fields  were  assumed  to  depend  only  on  depth, 
and  are  smooth  on  the  wavelength  scale.  The  linearized  inversion  is  designed  to  re¬ 
cover  the  wavelength  scale  fluctuations  (“reflectivity”)  in  bulk  modulus  and  density, 
treated  as  perturbations  of  the  background  fields. 

The  forward  and  adjoint  state  computations  were  based  on  a  finite  difference  simu¬ 
lator  for  viscoelastic  waves  developed  at  Rice  ([3]),  coupled  to  generic  optimization 
software.  Five  steps  of  the  conjugate  residual  algorithm  required  approximately  10 
CRAY  C90  CPU  hours,  reducing  the  normal  residual  in  each  experiment  to  roughly 
20%  of  its  initial  value.  Use  of  64  MW  of  memory  accommodated  20  snapshots,  which 
made  3  recomputations  using  Griewank’s  scheme  sufficient  to  compute  the  adjoint  in¬ 
volving  2272  timesteps.  The  total  of  floating  point  operations  was  increased  by  less 
than  100  percent  over  that  needed  by  a  “store  it  all”  algorithm. 

The  inverted  wavelength  scale  fluctuations  in  bulk  modulus  are  shown  in  Figure  1; 
the  analogous  density  result  is  very  similar.  Figure  2  and  3  show  cross  sections  of 
the  parameter  estimate  from  viscoacoustic  and  acoustic  inversion,  bulk  modulus  and 
density  reflectivity  respectively,  to  more  clearly  show  the  difference  between  the  results 
from  the  viscoacoustic  and  acoustic  inversion.  The  shallow  reflectors  are  the  same  for 
the  viscoacoustic  and  acoustic  inversion,  but  deeper  reflectors  are  underestimated  in 
the  purely  acoustic  inversion.  The  positions  of  the  reflectors  are  also  slightly  shifted 
due  to  dispersion  in  the  more  realistic  viscoacoustic  medium.  The  bulk  modulus  and 
density  reflectivity  are  for  the  most  part  completely  correlated.  Figure  4  shows  one 
of  the  few  places  in  the  estimates  where  they  are  not  correlated  (1800  m  depth  and 
840  m  offset).  There  is  a  zero  crossing  in  the  density  reflectivity,  with  hardly  any 
accompanying  wiggle  in  the  bulk  modulus  reflectivity.  This  anomaly  could  indicate 
a  significant  lithological  change.  The  result  is  very  preliminary,  since  the  normal 
equations  have  been  solved  rather  imprecisely. 
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Plane  1 
Trace 


Trace 
Plane  1 


Figure  1:  The  estimated  bulk  modulus  from  viscoacoustic  inversion.  Reflectivity  was 
not  estimated  above  1000  rn  depth,  roughly  corresponding  to  mute  of  the  data.  Depth 
scale  in  meters.  Aspect  ratio  approximately  1:1 


Figure  2:  A  depth  cross  section  of  the  estimated  bulk  modulus  reflectivity  at  750  m 
offset.  Solid  -  Estimate  from  viscoacoustic  inversion.  Dashed  -  Estimate  from  acoustic 
inversion. 
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Figure  3:  A  depth  cross  section  of  the  estimated  density  reflectivity  at  750  m  off¬ 
set.  Solid  -  Estimate  from  viscoacoustic  inversion.  Dashed  -  Estimate  from  acoustic 
inversion. 


Figure  4:  A  depth  cross  section  of  the  estimated  bulk  modulus  and  density  reflectivity 
at  840  m  offset.  Solid  -  Estimated  bulk  modulus  reflectivity.  Dashed  -  Estimated 
density  reflectivity.  Dash-dotted  encircled  area  contains  the  anomaly. 
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Conclusions  Griewank’s  algorithm  decreases  computational  cost  and  storage  re¬ 
quirements  dramatically  for  calculations  of  adjoints  of  non-conservative  systems.  The 
use  of  the  algorithm  enabled  a  viscoacoustic  inversion  of  a  real  data  set. 

Reflectivity  is  consistently  underestimated  and  incorrectly  located  if  anelastic  effects 
are  neglected.  Inversion  of  the  Mobil  AVO  data  set  revealed  a  tentative  reflectivity 
anomaly  below  a  prominent  reflector,  which  could  point  to  an  interesting  lithology. 

The  methods  developed  by  our  group  for  viscoacoustic  inversion  will  apply  as  well 
to  other  dissipated  transient  inverse  problems,  and  specifically  to  GPR  inversion. 
Interestingly,  radar  waves  penetrate  only  one  third  to  one  half  as  many  wavelengths 
into  the  subsurface,  even  under  ideal  conditions,  as  do  seismic  waves  at  exploration 
wavelengths.  Therefore  finite  difference  and  finite  element  based  inversion,  barely 
within  the  realm  of  feasibility  for  exploration  seismology,  should  be  a  viable  approach 
for  GPR  processing. 
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Abstract 

Recent  advances  in  inverse  scattering  algorithms,  coupled  with  the  availability  of  1-10 
GFLOP  computer  technology,  have  made  the  application  of  vector  electromagnetic 
inverse  scattering  methods  to  GPR  imaging  problems  practical.  Inverse  scattering 
methods  offer  the  promise  of  quantitative  reconstruction  of  the  3-D  variation  of  dielectric 
constant  and  conductivity  (and  possibly  the  magnetic  parameters)  in  the  earth. 

TechniScan  Inc.  (TSI)  has  developed  a  variety  of  algorithms  for  the  simulation  of  scalar 
and  vector  EM  scattering  problems  in  two  and  three  dimensions  [1-4].  These  algorithms 
include  frequency  domain  integral  equation  (FDIE),  finite-difference-time  domain 
(FDTD),  and  frequency  domain  recursion  (for  example,  one  way  parabolic  spatial 
marching)  methods.  Exact  Bessel  series  solutions  for  circular  cylindrical  (2-D)  objects 
buried  in  a  vertically  stratified  medium  illuminated  by  3-D  scalar  and  vector  fields  have 
also  been  developed  for  testing  the  accuracy  of  the  general  scattering  methods.  The  FDIE 
algorithms  are  particularly  advanced  and  have  been  successfully  applied  to  scalar  (TM 
polarization)  1-12.4  GHz  laboratory  GPR  data  for  2-D  cylindrical  targets  (Figures  1-5). 
Extension  of  the  algorithms  and  experimental  hardware  to  3-D  vector  GPR  is  the  major 
goal  of  our  present  research.  The  results  of  simulations  and  laboratory  experiments  will 
be  shown. 

A  critical  component  of  a  successful  GPR  imaging  approach  is  the  accurate  calibration  of 
the  antenna  fields  in  the  earth.  Inverse  scattering  methods  require  an  accurate  model  of 
the  magnitude  and  phase  of  the  applied  field  in  the  absence  of  the  targets  (incident  field) 
for  both  the  transmitting  and  receiving  antennas.  In  previous  2-D  work,  TSI  has 
developed  an  optimized  equivalent  source  (OES)  field  calibration  method.  In  this 
approach,  calibration  objects  of  known  properties,  both  in  air  and  buried  in  soil,  are 
scanned.  Equivalent  source  currents  postulated  on  a  surface  surrounding  the  antennas  (a 
line  segment  is  used  for  aperture  sources  in  2-D)  are  then  found  by  optimization  over  the 
calibration  data.  This  approach  has  been  verified  for  2-D  acoustic  and  electromagnetic 
data  collected  in  the  laboratory.  Extension  of  this  approach  to  3-D  vector  field  calibration 
is  an  important  goal  of  our  present  research.  Simulations  and  experimental  results  will  be 
shown. 

Another  problem  we  are  working  on  is  the  difficulty  involved  in  imaging  in  dispersive 
soils.  Sufficient  data  for  the  inversion  of  a  general  material  frequency  dependence  at 
each  image  pixel  is  not  available  in  reflection  mode.  TSI  has  developed  an  approach 
involving  the  application  of  the  Kramers-Kronig  causality  constraints  [5]  to  image 
frequency  dispersion  and  thus,  attenuation.  Details  of  this  approach  will  be  discussed  and 
preliminary  results  for  1-D  depth  profile  inversion  will  be  shown. 
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In  order  to  image  3-D  volumes  in  reasonable  time  using  O(GFLOP)  computers,  further 
advances  in  algorithm  speed  and  memory  efficiency  will  be  needed.  TSI  has  recently 
developed  a  3-D  parabolic  (one  way  wave  equation)  algorithm  for  transmission  mode 
(single  frequency)  imaging  of  low  contrast  acoustic  objects.  The  algorithm  has 

0(N3log(N))  complexity  for  one  transmitter  location  and  frequency.  Simulations  of  the 
female  breast  have  shown  sufficient  accuracy  and  speed  to  image  (100  X)3  volumes. 
Modifications  of  this  approach  have  been  shown  to  provide  accurate  backscattered  fields 
[6],  and  a  new  approach  [7]  uses  the  full  Helmholz  equation  and  thus,  removes  the  1-way 
restriction  completely.  For  time  domain  reflection  mode  imaging,  the  parabolic  approach 
requires  0(N4log(N))  computations  for  the'  evaluation  of  the  time  waveforms  for  one 
source  location.  We  are  examining  the  extension  of  this  approach  to  3-D  vector  EM 
reflection  mode  imaging.  Another  algorithm  option  to  examine  is  3-D  vector  FDTD 

which  has  the  advantage  of  providing  the  full  time  waveform  in  one  0(N4)  iteration.  We 
intend  to  compare  the  speed  and  accuracy  of  these  two  methods.  A  preliminary  analysis 
will  be  presented. 

Algorithm  and  experimental  research  plans  for  this  year  will  be  outlined. 
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Figure  1:  Experimental  geometry  of  the  laboratory  scale  GPR  scanner.  An  air-filled 
ABS  pipe  buried  in  sand  was  scanned  over  the  1-12.4  GHz  bandwidth  for  a  single 
transmitter/receiver  location  over  the  pipe.  This  is  sufficient  data  for  inversion  if  an 
angular  symmetry  constraint  of  the  target  is  applied.  A  1/2"  buried  alluminum  rod  was 
used  for  time  signal  calibration.  Cylindrical  wave  models  were  used  for  the  antenna  field 
spatial  variation  (in  air). 
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Figure  3:  Grayscale  of  the  true  Re(y  =  er-l)  (the  plastic  annulus  is  ignored  sine  its 
dielectric  constant  is  very  close  to  that  of  sand).  The  scale  is  0.2  (white)  to  -  0.7  (black). 


time  in  nanosec. 

Figure  2:  Comparison  of  the  lab  data  with  2-D  scalar  theory  for  the  experiment  shown  in 
Fig.  1.  Air-sand  vs.  pipe  amplitude  ratio  error  in  the  simulated  data  is  attributed  to 
inadequate  calibration  of  the  antenna  field  spatial  variation. 
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Figure  4:  Inverse  scattering  image  of  the  air-filled  pipe  from  lab  data. 
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Figure  5:  Centerline  profile  comparison  of  the  true  and  reconstructed  y  s  sr-l.  Accuracy 
of  the  image,  considering  the  air-sand/tube  signal  ratio  error,  is  attributed  to  the  angular 
symmetry  constraint,  coupled  with  the  accurate  arival  times.  Refinement  of  calibration 
methods  and  multiple  T/R  positions  should  obviate  the  need  for  this  constraint. 
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Abstract 

We  seek  a  thorough  understanding  of  the  factors  that  affect  the  recorded  response  of 
electromagnetic  and  seismic  wave  propagation  in  the  shallow  subsurface.  This  information 
will  help  us  to  evaluate  the  limitations  of  current  GPR  data  collection  and  processing 
methods  and  allow  us  to  investigate  ways  to  improve  image  resolution  and  target 
recognition.  A  typical  target  of  interest  in  this  research  would  be  a  contaminant  such  as  a 
DNAPL  (Dense  Non- Aqueous  Phase  Liquid). 

In  an  effort  to  understand  the  physical  contributions  to  ground  penetrating  radar  (GPR) 
response  of  contaminants  and  shallow  earth  media,  we  are  using  a  wave  propagation 
modeling  method  developed  in  the  frequency-wavenumber  domain.  Plane  waves  are 
propagated  through  a  layered  earth  described  by  variations  in  dielectric  permittivity, 
magnetic  permeability  and/or  electrical  conductivity.  The  ability  to  model  GPR  propagation 
in  realistic  earth  media  (dielectric  soils  with  conductive  losses)  will  help  us  investigate  ways 
to  improve  image  resolution  and  target  recognition  using  GPR  data.  We  also  hope  to 
understand  the  limitations  of  current  GPR  data  collection  methods  and  determine  the 
benefits  of  integrating  seismic  (acoustic  and/or  elastic)  data  with  the  radar  data. 

Maxwell’s  equations  of  electrodynamics  define  a  differential  system  that  can  describe  the 
propagation  of  plane  waves  in  lossy-dielectric  layered  media.  Exact  solutions  are  obtained 
for  the  electric  and  magnetic  field  intensities  by  solving  an  eigenvalue  problem  (ref.  1). 
The  local  characteristic  vectors  are  used  to  match  boundary  conditions  across  interfaces 
between  layers,  and  the  local  characteristic  values  are  the  permissible  vertical  wave 
numbers  within  the  medium.  Attenuation  due  to  non-zero  conductivity  is  included  directly 
through  the  complex  wavenumber.  This  is  not  as  easily  done  with  time-domain  methods 
(ref.  2).  Because  the  solutions  are  obtained  in  the  frequency  domain,  wave  dispersion 
effects  can  also  be  incorporated  through  the  permittivity,  permeability  and  conductivity 
functions.  A  fundamental  matrix  describes  propagation  through  the  layered  medium,  and 
the  reflection  and  transmission  response  of  the  medium  in  the  space-time  (recording) 
domain  is  determined  through  an  integration  of  the  fundamental  matrix  (refs.  3,4).  The 
matrix  method  described  allows  us  to  calculate  both  the  reflection  response  (reflectivity) 
and  the  transmission  response  (transmissivity)  of  the  layered  medium. 

When  modeling  monostatic  (vertically  incident)  radar  reflection  profiles,  we  can  include 
lateral  variations  in  the  layered  earth  model  along  the  profile  (a  different  model  for  each 
coincident  source-receiver  pair),  and  combine  signals  from  multiple  receivers  into  a  cross- 
section  image  of  the  subsurface  (refs.  5,6).  We  introduce  perturbations  to  the  background 
layered  medium  through  an  approximate  form  for  the  3-D  response  of  point  scatterers  (ref. 
7).  This  result  is  derived  from  the  calculation  of  an  analytic  form  for  the  convolution 
involving  the  scattered  wavefield  and  the  full  (background  plus  scattered)  wavefield. 
Building  up  the  response  of  (the  Green’s  function  for)  the  background  medium  one  layer  at 
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a  time,  the  response  of  any  point  scatterer  is  added  into  each  layer  calculation,  and  thus 
incorporated  directly  into  the  1-D  matrix  algorithm  described  above.  We  use  a 
superposition  of  the  scattered  wavefield  from  a  number  of  points  acting  as  Huygens’ 
sources  to  simulate  the  GPR  response  due  to  laterally  varying  media  (fig.  1). 

The  GPR  transmissivity  calculation  provides  an  indication  of  the  electrical  skin-depth  of  the 
layered  medium.  The  matrix  methods  above  allow  us  to  study  amplitude-as-a-function-of- 
depth  for  a  variety  of  GPR  source  characteristics  (frequencies,  angles  of  incidence)  and 
earth  parameters  (varying  conductivity,  dielectric  constant,  magnetic  permeability)  (figs. 
2,3).  Since  we  calculate  full  waveform  data,  we  can  also  study  the  phase  of  the  GPR 
signal  in  reflection  and  transmission  (refs.  6,8).  High  contrasts  in  conductivity  throughout 
the  earth  model  affect  the  reflection  and  transmission  coefficients,  and  this  produces  change 
in  both  the  amplitude  and  the  phase  of  the  traveling  waves  (fig.  4). 

Shallow  seismic  data  and  GPR  will  be  collected  along  coincident  profiles  during  a  field 
experiment  at  a  controlled  spill  test  site  at  Dover  AFB,  Delaware.  Seismic  signals  are 
sensitive  to  the  mechanical  properties  of  the  earth,  while  GPR  is  sensitive  to  the  subsurface 
electrical  properties.  GPR  operates  at  higher  frequencies  than  shallow  seismic,  but  does 
not  penetrate  as  far.  Seismic  surface  waves  may  inhibit  interpretation  of  the  shallowest 
body  wave  data,  however  this  region  is  well  sampled  by  the  GPR.  Through  analysis  of 
field  data  and  computer  modeling,  we  will  explore  ways  to  integrate  the  two  data  sets. 
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Figure  1.  Monostatic  GPR  reflection  profiles  over  two  simple  earth 
models,  (a)  GPR  at  1GHz  for  road-bed- with-void  model;  (b)  GPR  at  500 
MHz  for  buried-sand-channel-over-clay  model  with  DNAPL  simulated  by 
small  perturbation  in  the  propagation  constant. 
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Figure  2.  Transmission  amplitude  versus  depth  for  a  layered  earth  model 
where  conductivity  is  varied  in  alternating  layers  to  simulate  interbedded 
clays  and  sands.  At  zero  depth,  medium  grey  designates  unit  source 
amplitude  for  monochromatic  250  MHz  wavefield.  Lightening  of  the  grey 
shades  designates  reduction  of  amplitude  in  transmission  down  to  the  depth 
of  penetration  indicated  by  a  change  to  white  (1/e  level,  or  skin  depth). 


10  layers  alternating 


halfspace 


10 

7 

8 
9 


.001 

.01 

varied  from  .01  *  .06 

.01 

.01 


24 


depth  (in) 


Fixed  GPR  frequency  50  MHz 


conductivity  (mho/m) 


Interbedded  Layer  Earth  Mode! 

dielectric  conductivity 
constant  (mho/m) 


depth  (m) 


1.25 - - - 

halfspace  9  .01 


Figure  3.  Transmission  amplitude  as  ^function  of  depth  for 
monochromatic  50  MHz  source  wavefield.  Earth  model  and  grey-scale  as 
in  fig.  2. 
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Figure  4.  Reflection  and  transmission  GPR  records  after  varying  conductivity  in  interbedded  earth  model 
of  figures  2  and  3.  Pulse  radar  source  is  a  zero-phase  Ricker  wavelet.  Amplitude  and  phase  are  sensitive 
to  changes  in  conductivity  through  the  complex  permittivity. 
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Abstract 

Groundwater  flow  and  contaminant  transport  are  largely  governed  by  the  hydraulic 
conductivity  of  the  subsurface.  A  first  step  in  modeling  groundwater  flow  and/or 
contaminant  transport  at  a  site  is  to  obtain  a  3-D  model  of  the  hydraulic  conductivity  of  the 
subsurface.  Past  practice  has  often  involved  assigning  average  values  of  hydraulic 
conductivity  to  simple  layered  models  of  the  subsurface.  However  it  now  widely 
recognized  that  using  averaged  values  is  insufficient;  the  spacial  variability  in  conductivity, 
the  natural  heterogeneity  of  the  subsurface,  has  a  significant  effect  on  groundwater  flow 
and  is  a  critical  parameter  in  modeling  contaminant  transport.  There  are  two  key  issues 
associated  with  die  observed  variability  in  conductivity  that  are  the  focus  of  the  proposed 
research: 

1)  the  measured  hydraulic  conductivity  of  a  volume  of  the  subsurface  can  vary 
depending  upon  the  size  of  the  sampled  volume; 

2)  at  any  given  volume- scale  the  hydraulic  conductivity  exhibits  pronounced  spacial 

variability. 

Our  interest  in  this  proposal  is  in  addressing  both  of  these  issues  from  a  geophysical 
perspective.  Specifically,  we  plan  to  investigate  the  way  in  which  measurements  of 
dielectric  properties,  at  various  scales,  can  be  used  to  characterize  the  scale-dependence  and 
spacial  variability  of  hydraulic  properties. 

Methods  exist  that  allow  for  the  measurement  of  dielectric  properties,  at  high  resolution, 
and  across  a  wide  range  of  volume-scales.  Dielectric  properties  can  be  measured  in  the 
laboratory,  in  a  borehole  using  logging  tools,  in  a  borehole  or  across  a  2-D  exposure  using 
time  domain  reflectometry  (TDR),  at  the  surface  using  ground  penetrating  radar  (GPR), 
and  for  tomographic  imaging  using  borehole  antennas.  Given  the  wide  variety  of  methods 
available  for  in  situ  dielectric  measurement,  and  the  variation  in  sampled  volume-scale, 
there  is  tremendous  potential  for  using  dielectric  properties  to  characterize  the  hydraulic 
properties.  This  can  only  be  done  successfully  however  if  we  develop  a  clear 
understanding  of  the  link  between  dielectric  properties  and  hydraulic  properties,  across  a 
wide  range  of  volume-scales.  Of  particular  interest  in  this  proposal  is  the  use  of  ground 
penetrating  radar  as  a  means  of  obtaining  a  3-D  dielectric  "map"  of  the  subsurface.  The 
central  objective  of  the  proposed  research  is  to  determine  the  ways  in  which  the  measured 
dielectric  properties  can  then  be  related  to  unknown  hydraulic  properties. 

Our  first  step  has  been  to  define  a  terminology  for  referring  to  the  different  volume  scales. 
The  following  series  of  scales  is  adapted  from  those  used  in  modeling  fluid  flow  in  oil 
reservoirs1: 

microscopic  scale:  10's  to  100's  of  microns 
macroscopic  scale:  a  few  centimeters 
megascopic  scale:  10's  to  100's  of  meters 
gigascopic  scale:  kilometers 

The  microscopic  scale  is  the  size  of  a  few  pores.  At  this  scale  there  is  a  high  degree  of 
variability  in  the  properties  of  the  porous  medium,  and  the  continuum  equations  used  to 
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describe  flow  and  transport  are  not  considered  to  be  valid.  The  macroscopic  scale  is  the 
size  of  laboratory  samples  and  measurements,  on  the  order  of  two  to  five  centimeters.  It  is 
also  generally  assumed  that  this  is  the  scale  at  which  all  the  variability  seen  at  the  microscale 
is  averaged,  and  the  continuum  description  of  transport  becomes  valid.  Up-scaling  from 
the  macroscopic  scale,  we  reach  the  megascopic  scale.  This  scale  is  defined  as 
corresponding  to  the  size  of  hydrogeologic  units  such  as  sand  channels,  clay  lenses,  gravel 
beds.  The  largest  scale  is  the  gigascopic  scale,  which  is  defined  as  the  scale  of  the  aquifer. 

Our  research  approach  is  a  combination  of  laboratory  experiments,  theoretical  modeling  and 
field  experiments.  We  divide  the  research  problem  into  two  regimes:  those  length  scales 
(microscopic  and  macroscopic)  at  which  the  scale  of  dielectric  measurement  can  be  the 
same  as  the  required  conductivity  structure;  and  the  larger  volume-scales  (megascopic  and 
gigascopic)  at  which  the  scale  of  the  measured  dielectric  properties  is  much  smaller  than 
that  of  the  hydraulic  conductivity  structure. 

The  relationship  between  dielectric  and  hydraulic  properties  in  the  first  regime  is  addressed 
through  a  series  of  laboratory  studies  on  cylindrical  cells  of  standard  dielectrics  and 
unconsolidated  materials,  in  which  we  systematically  vary  the  material  properties  and  scales 
of  heterogeneity.  Measurements  are  made  of  dielectric  properties  and  hydraulic 
conductivity.  We  investigate,  at  the  laboratory  scale,  the  scale-dependence  of  each  property 
and  of  the  relationship  between  the  two.  These  results  will  have  implications  for  field  sites 
where  multi-frequency  GPR  can  be  used  to  obtain  dielectric  data  across  a  range  of  scales.  - 
Can  this  be  transformed  into  a  map  of  scale-dependent  hydraulic  conductivity  across  that 
same  range? 

In  our  laboratory  studies  we  are  working  with  a  silica  sand  and  a  standard  kaolinite  clay. 
The  sand  is  from  Lane  Mountain,  Washington,  and  is  essentially  pure  quartz.  The  clay  is 
from  the  Clay  Mineral  Repository  at  the  University  of  Missouri.  These  two  end-member 
materials  are  well  characterized  in  that  we  know  grain  densities,  grain  size  of  the  sand, 
particle  size  of  the  clay,  surface  areas,  and  cation  exchange  capacity  of  the  clay.  In  our 
previous  research  using  these  two  materials  we  have  developed  a  dielectric  model,  based  on 
the  complex  refractive  index  model  CREM^,  that  can  accurately  predict  that  variation  of 
dielectric  constant  with  porosity  and  clay  content.  We  start  by  defining  the  properties  of 
each  solid  component  in  its  "wetted"  state;  i.e.  we  determine  the  dielectric  constant  of  the 
solid  in  equilibrium  with  a  few  monolayers  of  adsorbed  water^.  Once  these  input 
parameters  are  known,  we  find  that  our  measured  data  can  be  well  modeled  using  the 
simple  volumetric  mixing  law,  CRIM: 


i 

where  Kt  is  the  dielectric  constant  of  the  total  measured  system  and  fi  and  Kj  are  the 
volume  fraction  and  dielectric  constant  of  the  ith  component 

From  our  experiments  with  the  sand/clay  mixtures  we  have  obtained  initial  results  that 
indicate  a  direct  link  between  dielectric  constant  and  permeability.  These  results  are  shown 
in  Figure  1  as  a  plot  of  measured  dielectric  constant  versus  permeability  for  a  suite  of 
sand/clay  mixtures.  The  dielectric  constant  was  measured  on  the  sand/clay  mixtures  at  a 
frequency  of  1  MHz  using  an  impedance  analyzer.  Permeability  was  measured  on  different 
sand/clay  mixtures,  but  with  the  same  composition,  using  a  falling  head  permeameter. 
There  is  very  good  agreement  between  our  predicted  relationship  and  that  observed  with 
our  data.  The  predicted  relationship  is  derived  using  CRIM  as  a  model  for  dielectric 
properties  and  the  Kozeny-Carmen  equation  as  a  model  for  permeability. 
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The  results  in  Figure  1  are  for  a  very  simple,  homogeneous  mixture  of  sand  and  clay.  Due 
to  the  constraints  of  the  measurement  techniques,  the  measurements  of  dielectric  constant 
and  permeability  were  made  on  different  samples,  with  different  sizes.  Our  present 
extension  of  this  research  is  to  conduct  experiments  on  heterogeneous  materials,  and  to 
assess  the  effect  of  sample  size.  We  have  recently  obtained  a  time  domain  reflectometry 
probe  which  can  be  used  with  the  same  samples  on  which  the  permeability  measurements 
are  made.  The  only  result  to  date  with  the  new  TDR  probe  is  a  measurement  of  the 
dielectric  constant  of  water;  not  a  new  result!  but  evidence  that  we  now  have  a  working 
and  accurate  probe.  We  are  presently  designing  the  samples  cells  for  the  joint 
permeability/dielectric  measurements. 

The  second  regime  of  scale  in  our  research  is  that  at  which  large-scale  hydrogeologic 
structures  are  assumed  to  dominate  the  flow.  This  is  the  scale  at  which  the  architecture  of 
the  aquifer  needs  to  be  characterized  in  terms  of  features  such  as  sand  channels,  gravel 
bars,  clay  lenses.  The  proposed  methodology  for  aquifer  characterization  at  this  scale 
involves  determining  the  sedimentological  environment  from  outcrop  and  core  and  then 
generating  a  sedimentological  facies  model  of  the  subsurface4.  Values  of  hydraulic 
conductivity  can  be  assigned  to  the  various  sedimentary  units,  thus  creating  a 
hydrogeologic  facies  model.  This  approach  has  great  potential  to  provide  realistic  models 
of  hydraulic  conductivity  variation,  as  it  is  based  on  years  of  research  by  sedimentologists 
who  have  studied  in  detail  the  spacial  variation  in  material  properties  within  different 
sedimentary  units.  A  limitation  to  this  method  is  the  need  to  adequately  characterize  the 
sedimentary  environment  in  order  to  generate  the  facies  model.  This  is  essentially 
impossible  to  do  in  areas  of  limited  outcrop  and  core  data. 

At  this  megascopic  scale,  there  is  enormous  potential  for  using  GPR  to  image  the  size, 
shape,  distribution  and  internal  structure  of  the  megascopic  hydrogeologic  units.  In  Figure 
2  is  shown  a  GPR  profile  from  an  area  east  of  Vancouver,  British  Columbia  where  we  are 
attempting  to  image  these  larger-scale  units.  As  shown  in  the  interpreted  section  there  is 
evidence  that  we  are  able  to  detect  cross-stratified  gravel  units,  clay  layers,  and  channels. 
The  use  of  GPR  for  identifying  such  larger-scale  units,  may  prove  to  be  an  extremely 
useful  way  of  obtaining  a  quantitative  description  of  the  mega-scale  heterogeneities 
controlling  flow  in  an  aquifer.  In  order  to  do  this  successfully  however,  there  need  to  be 
integrated  field  studies  that  allow  us  to  determine  the  electromagnetic  signature  of  various 
types  of  sedimentary  units;  i.e.  when  an  electromagnetic  wave  interacts  with  a  sand 
channel  in  a  stratified  sand/gravel  background  -  what  is  the  resulting  appearance  of  that 
feature  in  the  GPR  record  in  terms  of  the  velocity  structure,  and  the  amplitude  and 
continuity  of  reflectors? 

The  key  to  successful  imaging  and  thus  characterization  of  this  type  with  GPR,  is  the 
development  of  a  fundamental  understanding  of  which  hydrogeologic  features  are  "seen" 
with  GPR  and  why.  A  series  of  cliff  face  experiments  has  been  started  in  which  we  image 
features  visible  in  the  face  and  then  sample  in  detail  the  imaged  interface  to  determine  the 
dominant  parameters  responsible  for  the  successful  imaging.  Equally  important  will  be  the 
detailed  study  of  those  features  that  cannot  be  imaged  with  GPR.  We  conducted  our  first 
such  study  last  summer,  with  GPR  surveys  conducted  along  the  top  of  the  cliff  and 
samples  taken  for  measurement  of  porosity,  grain  size  and  moisture  content.  This  summer 
we  will  add  the  use  of  TDR  to  measure  in  situ  dielectric  constants.  We  are  proposing  to 
conduct  several  of  these  experiments,  imaging  specific  features  to  determine  which  can  be 
detected  with  GPR.  An  obvious  extension  of  these  experiments  is  the  use  of  forward 
modeling  to  attempt  to  simulate  the  recorded  GPR  data. 
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Having  completed  the  proposed  experiments  that  span  scales  from  centimeters  to  10's  of 
meters,  we  will  be  in  a  position  to  attempt  to  develop  a  conceptual  framework  for  the  scale- 
dependent  relationship  between  dielectric  and  hydraulic  properties.  At  this  point  in  this 
project  it  appears  that  the  two  defined  regimes  will  necessarily  correspond  to  two  separate 
regimes  for  theoretical  modeling  of  the  dielectric  and  hydraulic  behavior.  There  is  however 
the  intriguing  possibility  that  we  will  discover  a  means  of  scaling  the  physics,  so  that  we 
can  cross  the  boundary  from  macro-scopic  to  mega-scopic  in  a  way  that  integrates  a 
knowledge  of  geological  processes  to  predict  the  relationship  between  resulting  properties. 

In  conclusion,  the  objective  of  our  laboratory,  theoretical  and  field  study  is  to  investigate 
the  relationship  between  dielectric  properties  and  hydraulic  properties  across  a  wide  range 
of  volume-scales.  Our  approach  will  advance  the  understanding  of  the  scale-dependence  of 
these  properties,  and  more  importantly,  the  scale-dependence  of  the  relationship  between 
these  two  properties.  Ultimately  this  will  allow  us  to  more  effectively  use  remote 
measurement  of  dielectric  properties  as  a  means  of  characterizing  the  spacial  variability  in 
hydraulic  properties. 
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Figure  1:  Modeled  relationship  between  the  dielectric  constant  of  sand-clay  mixtures  and 
permeability  for  both  saturated  and  dry  materials.  The  measured  data  points  (observed) 
show  good  agreement  with  the  model  (predicted).  The  model  combines  a  complex 
refractive  index  model  of  the  dielectric  behavior  with  permeability  predictions  from  the 
Kozeny-Carmen  relationship. 
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Figure  2:  Example  of  GPR  data  collected  with  100  MHz  antennas.  Note  the  imaging  of  large-scale  hydrogeologic  features  as 
interpreted  in  the  lower  section. 
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Introduction  and  Background 

One  of  the  most  difficult  problems  in  environmental  remediation  is  location  of 
subsurface  organic  contaminants  [2,  18,  31,  32],  Most  common  organic  contaminants 
[8,  20,  22]  are  known  to  react  with  clay  minerals  [13,  Table  16,  and  references  therein]. 
However,  while  many  clay-organic  reaction  processes  are  known,  only  a  few  have  been 
the  subject  of  detailed  studies,  and  mostly  for  agricultural  purposes  [5,  7,  13,  21,  29]. 
Worse  yet,  some  of  these  reactions  alter  the  physical  properties  of  the  clays,  allowing 
contaminant  migration  through  supposedly  impermeable  clay  barriers  [4,  6,  10,  11,  12, 
34,  and  earlier  references  reviewed  in  13,  p.  516-525). 

Although  the  induced  polarization  (IP)  process  was  observed  by  Conrad  Schlumberger 
as  early  as  1912,  it  was  not  until  the  1950's  that  mining  companies  used  the  induced 
polarization  (or  complex  resistivity)  method  to  detect  economic  minerals  (see  reviews  in 
3,  28]).  Complex  resistivity  does  this  by  detecting  the  electrochemical  processes  that 
occur  between  minerals  and  water  (such  as  oxidation-reduction  corrosion  of  sulfides, 
see  review  in  [15,  17]).  Petroleum  companies  have  exploited  IP  effects  caused  by 
hydrocarbon  seepage  that  results  in  iron  sulfide  geochemical  alteration  halos  [26].  More 
recently,  complex  resistivity  has  been  exploited  to  observe  cation  exchange  at  clay- 
water  surfaces  and  to  directly  measure  the  reactions  between  organics  and  clay  minerals 
(see  reviews  in  [15, 18]). 

In  1984,  the  U.  S.  Environmental  Protection  Agency  Environmental  Monitoring 
Systems  Laboratory  and  the  U.  S.  Geological  Survey  Water  Resources  Division  were 
investigating  organic  contamination  problems  in  Nevada  and  Louisiana,  and  requested 
geophysical  assistance.  At  that  time,  the  detection  of  organics  with  geophysics  was  not 
at  all  understood  nor  expected.  Core  samples  of  representative  uncontaminated  and 
organic-contaminated  materials  from  the  sites  were  measured  and  a  large  IP  difference 
was  observed  in  the  laboratory  between  the  uncontaminated  and  contaminated  materials 
[16].  Although  the  cause  of  the  IP  response  was  not  known  nor  understood,  it  was  large 
enough  to  be  exploited  in  the  noninvasive  mapping  of  organics  in  the  subsurface[16].  It 
also  showed  that  at  the  Nevada  site,  reaction  of  the  organic  with  a  impermeable  clay 


*  1500  Illinois  St.,  Golden,  CO  80401-1887,  (303)  273-3458,  -3478  fax, 
golhoeft@mines.colorado.edu;  formerly  with  the  U.  S.  Geological  Survey,  Denver,  CO 


33 


barrier  supporting  a  perched  water  table  had  caused  the  contaminants  to  penetrate 
through  more  than  a  meter  of  100%  montmorillonite  clay  [16].  The  clay  layer  was  of 
the  swelling  type  that  would  squeeze  a  drill  hole  shut  in  a  few  days,  so  the  penetration 
by  the  organics  was  a  suiprise.  Unfortunately  at  the  Nevada  site,  the  contaminant  plume 
was  formed  from  many  disparate  sources  and  a  large  variety  of  organic  chemicals, 
which  have  reacted  amongst  themselves  to  create  new,  not-before-known  organics,  and 
it  has  not  been  possible  to  determine  which  ones  are  specifically  causing  the  clay- 
organic  response. 

The  Louisiana  site  was  a  much  simpler  mix  of  contaminants  resulting  from  the  disposal 
of  oil  field  brines  and  hydrocarbons  in  a  landfill  [16].  To  understand  and  better  exploit 
the  field  observations,  a  systematic  laboratory  investigation  was  begun  to  identify  the 
organics  reacting  with  the  clays  [23]  and  to  identify  the  mechanisms  and  reaction 
processes.  Toluene  was  identified  as  the  reactive  organic  at  the  Louisiana  site,  allowing 
duplication  of  the  field  observations  and  core  studies  with  synthetic  mixtures  of  toluene 
and  montmorillonite  [23]  using  clay  mineral  standards  [33].  Further  laboratory 
investigations  [9,  19]  found  the  reaction  mechanism  to  be  the  polymerization  of  toluene 
to  bibenzyl  in  the  Lewis  Acid  surface  environment  of  montmorillonite  clay.  Identifying 
and  understanding  the  mechanism  of  the  clay-organic  (montmorillonite-toluene) 
reaction  (polymerization  to  bibenzyl)  allowed  calibrating  the  geophysical  field 
measurements  and  turning  anomaly  mapping  into  quantitative  toluene  concentration 
mapping  (Olhoeft,  in  press). 

In  addition  to  mapping  organic  contamination  by  their  reactions  with  clays,  complex 
resistivity  should  also  be  able  to  monitor  changes  in  the  properties  of  the  clays  caused 
by  reaction  with  the  organics.  As  mentioned  and  referenced  above,  such  processes  can 
have  adverse  impact  on  contaminant  migration  by  causing  the  breach  of  otherwise 
impermeable  clay  barriers.  Thus,  complex  resistivity  could  be  a  useful  monitoring  tool 
to  periodically  test  the  integrity  of  landfill  clay  mineral  liners  or  natural  clay  mineral 
impermeable  barriers. 

Finally,  the  complex  resistivity  monitoring  of  clay-organic  processes  could  also  be 
useful  in  monitoring  unanticipated  and  undesirable  changes  during  a  site  remediation. 
In  particular,  it  should  be  an  excellent  tool  to  monitor  remediation  using  electrokinetic 
effects  [1,  24,  25]. 

Research  Effort 

Synthetic  mixtures  of  clay  minerals  and  selected  organic  contaminants  will  be  measured 
in  the  laboratory  to  determine  which  combinations  produce  observable  complex 
resistivity  responses  as  a  function  of  frequency.  Those  with  observable  responses  will 
be  further  tested  as  a  function  of  clay  mineral  type,  water  chemistry  and  organic 
concentration  to  determine  the  limits  of  the  complex  resistivity  response.  Those  with 
useful  response  limits  will  be  investigated  to  determine  the  mechanisms  and  processes 
of  clay-organic  interaction  using  standard  geochemical  laboratory  procedures. 
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The  details  of  this  investigation  will  have  to  evolve  as  the  investigation  proceeds  as  it  is 
basic  research  into  the  unknown.  In  the  beginning,  mixtures  of  single  components  of 
most  common  organic  contaminants  [13]  will  be  mixed  with  standard  and  well 
characterized  clay  minerals  [33]  and  the  complex  resistivity  response  will  be  measured. 
Cross-discipline  expertise  in  Geophysics,  Geology,  Chemistry  and  Geochemistry, 
Petroleum  Engineering  and  in  Environmental  Sciences  will  be  required  to  follow 
whatever  path  the  clay-organic  investigation  requires.  In  the  case  of  the 
montmorillonite-toluene  reaction  [19],  diffuse  infrared  reflectance  spectroscopy  proved 
the  most  useful  tool  in  unravelling  the  process.  In  other  cases,  different  tools  may  prove 
more  useful  and  a  comprehensive  suite  of  investigatory  techniques  is  required  and 
available  through  interdisplinary  cooperation. 

Natural  soil  samples  from  the  Air  Force  testbed  at  Dover,  Maryland  will  also  be  run 
through  the  laboratory  investigations  to  determine  the  possible  clay-organic  interactions 
at  the  site  for  the  organic  contaminants  to  be  tested  at  that  site.  If  the  site  materials 
show  promising  complex  resistivity  response,  noninvasive  field  complex  resistivity 
measurements  will  be  performed  to  map  the  geochemical  reactive  heterogeneity  of  the 
site  and  provide  a  baseline  for  future  testing. 
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Abstract 

Imaging  methods  are  based  mostly  on  small  variation  representations  of  the  unknown 
parameters  of  the  medium  (the  Born  representation).  This  is  done  extensively  in  cross¬ 
hole  electromagnetic  imaging  of  the  earth  (D.  L.  Alumbaugh,  Berkeley  Thesis  1993)  and  in 
impedance  tomography  in  general.  In  traveltime  tomography  Berryman  (Inverse  Problems, 
6,  21-42,  1990)  introduced  constrained  reconstruction  algorithms  based  on  Fermat’s  varia¬ 
tional  principle  that  work  even  with  high-contrast  media.  This  idea  extends  to  impedance 
tomography  whenever  variational  principles  are  available  as  in  the  static  or  very  low  fre¬ 
quency  cases. 

In  this  lecture  I  will  discuss  the  use  of  variational  principles  to  analyze  direct  and  inverse 
problems  with  high  contrast  in  the  low  frequency  or  quasistatic  regime  which  is  most 
useful  for  underground  imaging.  I  will  present  extensive  numerical  results  that  show  how 
high-contrast  difficulties  can  be  approached  and  sometimes  overcome. 

This  is  joint  work  with  J.  Berryman  and  L.  Borcea. 

Acknowledgement /Disclaimer 

This  work  is  sponsored  by  the  Air  Force  Office  of  Scientific  Research  under  Grant  F49620- 
94-1-0436.  The  views  and  conclusions  contained  herein  are  those  of  the  authors  and  should 
not  be  interpreted  as  necessarily  representing  the  official  policies  or  endorsements,  either 
expressed  or  implied,  of  the  Air  Force  Office  of  Scientific  Research  or  the  US  Government. 


38 


MODELING  AND  INVERSION  OF  SHALLOW  SEISMIC 
DATA  INCLUDING  NONGEOMETRICAL  WAVES 

John  Scales  and  Ilya  Tsvankin 
Center  for  Wave  Phenomena 
Colorado  School  of  Mines,  Golden 


Abstract 

Our  research  is  focused  on  advanced  seismic  inversion  and  imaging  methods  based  on  a 
more  comprehensive  treatment  of  wave  propagation  than  the  classical  ray  theory.  A  ma¬ 
jor  component  of  the  project  is  investigation  of  the  so-called  “nongeometrical”  waves  in 
realistic  models  of  the  subsurface  (the  term  “nongeometrical”  refers  to  waves  that  cannot 
be  accounted  for  by  the  geometrical-seismics  approximation  conventionally  used  in  seis¬ 
mology).  Nongeometrical  waves  can  supply  information  about  some  parameters  of  the 
medium  that  are  poorly  constrained  by  conventional  waves,  such  as  the  elastic  properties 
of  waveguides.  The  information  contained  in  full  waveforms  makes  it  possible  to  devise 
enhanced  methods  of  imaging  of  complicated  geologic  structures  and  different  types  of 
inhomogeneities. 

To  take  full  advantage  of  this  comprehensive  treatment  of  seismic  wavefields,  we  are 
developing  inversion  methods  capable  of  producing  subsurface  images  with  higher  accuracy 
and  resolution  than  conventional  algorithms. 


Fig.  1. 

We  have  adopted  a  multi-faceted  approach  to  full- waveform  inversion  including:  (1)  in¬ 
vestigation  of  novel  objective  functions  which  reduce  the  extreme  non-convexity  of  full- 
waveform  inverse  problems;  (2)  a  new  statistical  approach  to  incorporation  of  realistic  prior 
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information;  (3)  effective  combination  of  global  optimization  methods  with  rapid,  gradient- 
descent  algorithms.  We  have  begun  developing  an  extensive  modeling  capability  (based  on 
finite-difference  methods)  on  the  massively  parallel  Connection  Machine  (CM-5).  We  are 
also  working  on  a  library  of  nonlinear  local  and  global  optimization  routines,  including  a 
variety  of  Monte  Carlo  methods,  nonlinearly  constrained  downhill  simplex,  as  well  as  other 
direct  search  methods. 

The  advantages  of  our  approach  should  be  especially  significant  in  shallow  surveys,  which 
are  strongly  affected  by  near-field  phenomena  and,  consequently,  nongeometrical  waves. 
The  level  of  resolution  to  be  achieved  with  our  new  technologies  will  depend  on  the  probing 
wavelength  and  the  scattering/attenuative  properties  of  the  subsurface.  While  it  is  difficult 
to  give  a  quantitative  estimate  of  the  resolution  before  the  work  is  complete,  we  expect 
a  significant  improvement  over  conventional  methods  because  our  algorithms  will  treat 
the  whole  wavefield  (including  nongeometrical,  refracted  waves  etc.)  rather  than  just  the 
first  arrivals  and  make  no  smoothing  assumptions  about  the  medium.  We  expect  that  our 
methods  will  be  helpful  in  different  aspects  of  site  characterization  and  monitoring,  such  as 
subsurface  contaminant  location,  groundwater  studies,  and  other  environmental  problems. 

Our  project  is  currently  at  its  initial  stage.  In  the  near  future  we  will  be  working  on 
full-waveform  modeling  and  inversion  of  cross-hole  seismic  data  with  the  goal  of  obtain¬ 
ing  high-resolution  seismic  images.  In  addition  to  synthetic  models,  we  will  apply  our 
algorithms  to  a  multicomponent  cross-hole  data  set  recorded  by  Conoco  at  its  Newkirk, 
Oklahoma,  test  site.  This  site  is  characterized  by  a  vertically  stratified  geology  show¬ 
ing  large  velocity  contrasts.  This  combination  results  in  complicated  raypaths,  guided 
wave  phenomena,  and  possible  generation  of  nongeometrical  waves.  Figure  2  is  a  sample 
three-component,  common-receiver  gather  showing  significant  energy  in  late  arrivals.  The 
horizontal  components  contain  intensive  low-frequency  guided  waves  marked  by  B.  It  is 
likely  that  other  nongeometrical  waves  contribute  to  the  wavefield  between  the  first  arrivals 
and  guided  modes  (A  in  Figure  2).  Obviously,  these  are  only  preliminary  observations 
which  should  be  supported  by  a  more  detailed  investigation. 

Although  we  are  not  planning  on  studying  other  geophysical  techniques  within  the  frame¬ 
work  of  the  project,  the  inversion  methods  we  are  working  on  will  be  generally  applicable 
in  a  wide  variety  of  imaging  technologies.  For  instance,  our  inversion  algorithms  can  be 
used  to  enhance  radar  detection  methods  by  properly  modifying  the  objective  functions. 
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Fig.  2.  Three-component  common-receiver  gathers  for  a  single  source  polarization.  The  left 
column  is  the  vertical  displacement  component,  the  next  two  columns  represent  horizontal 
components.  A  and  B  are  explained  in  the  text.  Data  courtesy  of  Conoco. 
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Introduction 

High-frequency  (100  MHz  to  1  GHz)  ground-penetrating  radar  methods  have  shown 
excellent  resolving  capabilities  at  very  shallow  depths.  The  depth  of  penetration  of  radar 
energy,  however,  varies  widely  depending  on  the  soil  electrical  properties,  and  is  generally 
too  shallow  for  many  environmental  site  characterization  surveys.  A  promising  technique  for 
imaging  shallow  buried  targets  is  electromagnetic  induction.  In  order  to  have  reasonable 
penetration  into  the  subsurface,  we  must  use  relatively  low  frequencies,  e.g.,  a  few  kHz  to  a 
few  tens  of  MHz.  The  equation  that  applies  to  electromagnetic  induction,  over  most  of  the 
frequency  range  of  interest,  is  actually  the  diffusion  equation,  rather  than  the  wave  equation. 
Objects  are  detected  in  the  subsurface  via  secondary  fields  induced  in  the  buried  objects  rather 
than  propagating  energy  being  reflected  off  the  objects. 

A  number  of  commercial  electromagnetic  instruments  are  currently  used  for  shallow 
subsurface  electromagnetic  induction  surveys.  These  systems  are  effective  for  determining 
depths  to  layers  (for  example,  depth  to  water  table  or  depth  to  bedrock)  and  for  detection  of 
buried  targets.  They  have  had  limited  success  for  imaging  of  buried  targets  in  the  subsurface. 

The  LASI  High-Resolution  EM  Systems 

The  University  of  Arizona,  Laboratory  for  Advanced  Subsurface  Imaging  (LASI)  has 
developed  research  electromagnetic  (EM)  sounding  systems  that  provide  a  significant  advance 
in  the  state-of-the-art  for  shallow  subsurface  imaging  (Sternberg  et  al.,  1991;  Sternberg  and 
Poulton,  1994).  These  systems  record  the  ellipticity  of  the  received  magnetic  field  signal. 
The  received  magnetic  field  from  a  remote  transmitter  traces  an  ellipse  as  a  function  of  time. 
The  ellipticity  is  defined  as  the  ratio  of  the  minor  axis  of  the  ellipse  to  the  major  axis. 
Ellipticity  is  a  particularly  sensitive  parameter  for  mapping  depths  to  buried  layers  in  the 
subsurface  (Hoversten,  1981). 

Figure  1  shows  a  photograph  of  the  LASI  high-resolution  EM  transmitter.  The  electronics 
are  mounted  on  a  6-wheel  drive  all  terrain  vehicle  (ATV).  The  transmitting  coil  is  mounted 
on  a  boom  extended  in  front  of  the  ATV.  Figure  2  shows  a  photograph  of  the  receiver  ATV. 
The  receiver  coil  is  mounted  on  the  boom  at  the  rear  of  the  ATV.  This  EM  system  measures 
ellipticity  over  frequency  ranges  of  1  kHz  to  1  MHz  or  32  kHz  to  32  MHz.  An  earlier  version 
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of  this  system  recorded  over  a  frequency  range  of  30  Hz  to  30  kHz.  These  systems 
incorporate  a  number  of  unique  features,  including:  (1)  high-accuracy  simultaneous 
calibration,^)  automated  data  acquisition,  (3)  high  signal/noise  ratio,  and  (4)  rapid  surveying 
with  very  dense  measurements.  The  LAST  systems  represent  a  significant  advance  in  the 
state-of-the-art  for  shallow  EM  sounding. 

These  EM  sounding  systems  form  the  foundation  for  the  new  methods  that  are  being 
developed  on  this  project. 

Testing  and  Demonstration 

Initial  testing  and  demonstration  of  the  imaging  techniques  developed  in  this  project  will  be 
carried  out  at  a  physical  modeling  facility  that  we  recently  built  at  our  Avra  Valley  Test  Site. 
The  upper  end  of  the  frequency  range  in  the  EM  imaging  systems  will  include  effects  from 
both  conduction  currents  and  displacement  currents.  It  is  difficult  to  obtain  reliable  numerical 
modeling  calculations  of  theoretical  responses  to  complex  targets  in  this  frequency  range.  We 
will  therefore  use  full-size  physical  models. 

The  modeling  tank  is  located  at  our  test  site  in  Avra  Valley,  Arizona,  west  of  the  University 
of  Arizona  campus.  The  tank  is  20  m  long  by  3  m  deep  by  6  m  wide  (see  Figure  3).  The 
transmitter-receiver  array  will  be  kept  stationary  in  order  to  avoid  variations  in  response  due 
to  background  effects.  Various  targets  will  be  moved  in  the  tank  along  a  profile  line  under 
the  array.  Repeated  measurements  will  be  made  with  the  targets  at  different  depths, 
orientations,  and  with  different  types  of  targets.  The  background  conductivity  and  dielectric 
constant  of  the  fluid  in  the  tank  will  be  varied  by  using  mixtures  of  different  fluids.  Only  non- 
hazardous  fluids  will  be  used.  This  physical  modeling  facility  provides  a  unique  capability  for 
testing  and  demonstration  of  the  new  EM  imaging  techniques. 

Final  testing  of  the  EM  imaging  techniques  developed  in  this  project  will  be  carried  out  at  a 
lined  basin  which  is  also  located  at  our  Avra  Valley  Test  Site.  The  basin  was  designed  to  test 
the  ability  of  the  high-resolution  EM  system  to  image  the  flow  of  fluid  in  the  subsurface.  It 
is  approximately  5  m  deep,  and  30  m  x  30  m  square.  Heavy-duty  high-density  polyethylene 
(HDPE)  plastic  was  used  to  line  the  bottom  and  sides  of  the  excavated  basin.  Gravel  was 
then  placed  in  the  bottom  of  the  basin  along  with  drain  pipes.  Figure  4  shows  a  photograph 
of  the  liner  and  the  completed  pit  after  the  gravel  layer  had  been  placed  in  the  bottom  and 
before  soil  had  been  filled  in.  The  lined  basin  provides  a  closed  system  for  injection 
experiments.  It  is  a  unique  facility  for  monitoring  the  flow  of  fluids  which  provides  minimal 
interference  from  the  bottom  and  sides  of  the  pit. 

EM  Sounding  with  the  LASI  Systems 

Figures  5  and  6  show  examples  of  EM  soundings  with  the  LASI  high-resolution  EM 
Sounding  system  (Sternberg,  1993).  This  example  used  the  lined  basin  shown  in  Figure  4. 

During  the  summer  of  1992,  24,170  liters  of  fluid  were  injected  along  a  1  m  x  25  m  strip 
through  the  middle  of  the  lined  basin.  High-resolution  subsurface-imaging  ellipticity  surveys 
were  made  over  the  injection  region  prior  to  and  during  the  injection.  Figure  5  shows  an 
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image  of  the  electrical  resistivity  of  the  ground  immediately  after  the  injection  started.  Figure 
6  shows  the  electrical  resistivity  of  the  ground  1 7  days  after  the  start  of  the  injection.  These 
images  produced  by  the  EM  system  agree  very  well  with  cross-sections  of  ground  resistivity 
as  determined  by  well-log  measurements  that  were  made  in  25  boreholes  which  span  the 
injection  region.  The  EM  imaging  system  was  very  successful  in  imaging  the  location  of  the 
fluid  during  the  entire  course  of  the  injection.  Note,  however,  that  this  fluid  plume  is 
relatively  large  and  extended  effectively  infinite  in  and  out  of  the  section. 

Current  Research 

Ellipticity  measurements  contain  information  about  the  average  properties  of  the  earth 
between  the  transmitter  and  receiver.  Thus,  ellipticity  is  best  suited  for  measuring  the 
properties  of  a  layered  earth.  Tests  have  shown  ellipticity  measurements  to  be  relatively 
insensitive  to  small  three-dimensional  targets.  In  order  to  detect  relatively  small  3D  targets, 
we  have  redirected  our  current  research  toward  improving  the  imaging  capabilities  of  these 
EM  systems.  Specifically,  we  are  studying  arrays  of  magnetic  field  sensors,  including:  (1) 
null-coupling  arrays  for  greater  sensitivity,  (2)  focusing  arrays  for  greater  resolution,  (3) 
three-dimensional  arrays  for  improved  target  identification,  and  (4)  large,  dense  arrays  for 
faster,  lower-cost  acquisition. 

The  initial  results  with  prototype  arrays  are  very  encouraging.  Greatly  improved  sensitivity 
and  resolution  appear  to  be  obtainable.  Some  fundamental  problems  remain,  including:  (1) 
reducing  the  capacitive  coupling  between  the  transmitter  and  receiver,  and  (2)  determining 
optimum  array  sizes  and  configurations. 

Acknowledgment/Disclaimer 

This  work  was  sponsored  (in  part)  by  the  Air  Force  Office  of  Scientific  Research,  USAF, 
under  grant/contract  number  F49620-95- 1-0004.  The  views  and  conclusions  contained 
herein  are  those  of  the  authors  and  should  not  be  interpreted  as  necessarily  representing  the 
official  policies  or  endorsements,  either  expressed  or  implied,  of  the  Air  Force  Office  of 
Scientific  Research  or  the  U.S.  Government. 

References 

1 .  G.  M.  Hoversten,  Comparison  of  Time  and  Frequency  Domain  EM  Sounding  Techniques, 
Ph.D.  Thesis,  University  of  California,  Berkeley,  1981. 

2.  B.  K.  Sternberg,  S.  J.  Thomas,  N.  H.  Bak  and  M.  M.  Poulton,  High-Resolution 
Electromagnetic  Imaging  of  Subsurface  Contaminant  Plumes,  EPRI EN-7519,  1991. 

3 .  B.  K.  Sternberg,  Construction  of  a  Lined  Basin  for  Tests  of  a  High-Resolution  Subsurface 
Imaging  Ellipticity  System,  EPRI  TR- 103 462,  1993. 

4.  B.  K.  Sternberg  and  M.  M.  Poulton,  High-Resolution  Subsurface  Imaging  and  Neural 
Network  Recognition,  Proc.  SAGEEP  '94,  Boston,  1994,  pp.  847-855. 


44 


APPARENT  DEPTH  (METERS)  APPARENT  DEPTH  (METERS) 


A  PDE  BASED  APPROACH  TO  DISPERSIVE  INVERSE  PROBLEMS 


Yun  Wang 

Mathematical  Products  Division 
Armstrong  Laboratory,  OES 
Brooks  AFB,  TX  78235 


Abstract 

A  time  domain  approach  for  study  of  electromagnetic  inverse  problems  is  presented 
in  this  note.  Maxwell  equations  coupled  with  a  generalized  dispersion  model,  an 
electric  polarization  equation  governed  by  an  n-th  order  differential  equation  are 
considered.  Well  posedness  is  presented  for  a  one-dimensional  dispersive  medium 
case.  Parameters  representing  the  electromagnetic  property  of  a  medium  are  the 
static  permittivity,  relaxation  time,  natural  frequency,  static  conductivity,  etc..  The 
number  of  parameters  involved  depend  on  the  choice  of  order  of  derivatives  in  the 
polarization  equation. 

In  environmental,  clinical  medicine  and  many  other  areas,  the  microwave  images  of 
tissue  structures  and  soils  play  very  important  roles.  The  microwave  images  are  useful 
in  detection/remediation  of  underground  toxic  wastes,  and  detection/enhanced  treat¬ 
ment  of  abnormality  of  human  organs  and  tissue.  The  electromagnetic  properties  of 
a  medium  are  characterized  by  its  electric  polarization  mechanisms  and  its  static 
conductivity.  In  this  note  we  focus  on  the  development  of  partial  differential  equa¬ 
tions  (Maxwell  equations)  based  identification  techniques  for  physical  and  biological 
distributed  parameter  systems,  with  those  for  living  tissue  being  a  special  case.  We 
attempt  to  estimate  the  conductivity  and  parameters  which  characterize  the  polariza¬ 
tion  of  living  tissue  based  on  the  incident  and  scattered  electromagnetic  signals.  We 
point  out  that  in  literature  the  complex  permittivity  and  conductivity  which  depend 
on  frequency  of  the  emitting  signal  are  commonly  used  as  electromagnetic  proper¬ 
ties.  The  functions  of  complex  permittivity  and  conductivity  are  determined  by  a 
polarization  mechanism  which  is  assumed  to  be  governed  by  an  n-th  order  ordinary 
differential  equation  (see  [11,  16]  for  the  case  n  =  2).  Since  our  approach  is  in  the 
time  domain,  we  focus  directly  on  the  polarization  equations. 

The  electromagnetic  inverse  problem  has  been  studied  for  at  least  two  decades.  A 
good  survey  of  existing  methods  is  given  by  Albanese  et  al.  in  [1].  However,  there  is 
very  little  literature  on  the  study  of  electromagnetic  inverse  problems  in  the  time  do¬ 
main  employing  variational  formulations.  The  variational  formulation  approach  has 
been  successfully  applied  to  damped  hyperbolic  systems  in  [5,  6,  7]  and  hybrid  systems 
in  [4] .  A  similar  approach  applied  to  Maxwell  equations  can  be  found  with  focus  on 
well  posedness  and  control  problems.  For  example,  Duvaut  and  Lions  [9]  proved  ex¬ 
istence  and  uniqueness  of  Maxwell  equations  for  a  three-dimensional  inhomogeneous 
medium  with  superconductive  boundary.  The  medium  is  stable  (nonzero  static  con¬ 
ductivity)  and  polarization  is  assumed  to  be  proportional  to  the  electric  field.  Other 
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efforts  focused  on  controllability  and  stabilization  in  the  context  of  semigroups  and 
variational  formulations  can  be  found  in  [3,  12,  13,  14,  15]. 


We  consider  different  dispersion  models  (i.e.  models  for  polarization)  in  the  context 
of  Maxwell  equations  in  a  theoretical  study  of  the  type  and  rate  of  attenuation  of 
signals  associated  with  a  given  dispersion  model  for  living  tissue. 


For  a  three-dimensional  inhomogeneous  medium,  the  equation  governing  electromag¬ 
netic  phenomena  are  Maxwell  equations 

VD  =  0 

VB  =  0 

„  „  dB 

VxE  =  ~~§t 

VxH  =  ^2  +  Ji+Jo 
at 

and  constitutive  relations 

(i) 

Jc  =  <7  E 

(2) 

D  —  6o  E  +  P 

(3) 

B  =  no  H  +  M, 

(4) 

where  E  is  the  electric  field  intensity,  D  is  the  electric  flux  density,  H  is  the  magnetic 
field  intensity,  B  is  the  magnetic  flux  intensity,  J;  is  the  source  electric  current  den¬ 
sity,  P  is  the  electric  polarization,  M  is  the  magnetic  polarization,  a  is  the  position 
dependent  static  conductivity,  eo  and  Ho  are  the  permittivity  and  permeability  in 
vacuum  respectively.  The  bold  faced  characters  are  vectors  in  Cartesian  coordinates. 

In  preliminary  research,  the  propagation  of  a  plane  wave  is  investigated.  Assuming 
the  plane  wave  is  to  be  uniform  in  planes  parallel  to  the  x-y  plane  and  propagates  in 
z  direction,  then  the  electric  and  magnetic  field  intensity  are  reduced  to 

E  =  E(t ,  z )  in  x  direction 

H  =  H (t,  z)  in  y  direction. 


The  scalar  fields  E  and  H  will  be  adopted  from  now  on  with  understanding  the  E  is 
polarized  in  x  direction  and  H  is  polarized  in  y  direction.  Furthermore,  we  assume 
that  magnetic  polarization  is  zero,  i ,e.,M(t,z)  =  0  which  is  a  good  approximation  to 
human  tissue  [2].  With  the  above  assumptions,  Maxwell  equations  coupled  with  the 
constitutive  relations  yield  a  second  order  equation 


d2E  dE  d2P  d2E  dJs 

'“’£°  W  +  ''o£°  Ti  + W  -  _  W’ 


t>  0,  0<z<l.  (5) 
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While  tissue  is  thought  to  have  magnetic  properties  of  vacuum,  dispersion  of  electro¬ 
magnetic  signals  in  heterogeneous  media  is  a  complex  phenomenon  which  is  usually 
accounted  for  in  the  polarization  vector  P.  In  the  case  of  dispersive  media  (here, 

we  adopt  the  definition  by  Stratton  [16]  which  states  that  a  medium  is  said  to  be 

dispersive  if  the  phase  velocity  in  the  medium  is  a  function  of  frequency),  we  consider 
a  general  representation  for  the  electric  polarization  described  by  n-th  order  ordinary 
differential  equation 

dnp  «  Qn-ip 

aF  +  ^t  ’W=*=a°E'  (6) 

where  the  coefficients  are  constants.  This  representation  of  polarization  takes 

into  account  the  molecular  constitution  of  matter  and  treats  the  molecules  as  dynam¬ 
ical  systems  possessing  natural  frequency.  For  the  case  of  n  =  0, 1,2  the  polarization 
mechanisms  have  been  studied  extensively  [11,  16].  When  n  =  0,  the  nonzero  con¬ 
stant  do  is  identified  as  the  product  of  electric  susceptibility  xe  and  eo,  hence  the 
electric  flux  density  (3)  becomes 

D  =  (1  +  xe)coE  =  eE , 


in  which  t  is  called  the  static  permittivity  of  the  medium.  Many  media  (such  as 
metals)  exhibit  this  property. 


When  n  =  1,  we  can  identify  (6)  with  the  Debye  model  (which  is  considered  to  the 
polarization  mechanism  of  water  [1]) 

dP 

~E7  +  o\P  =  ao  E 


which  is  commonly  expressed  as 


dP 


E 


where  es  is  the  static  relative  permittivity  and  is  the  limiting  permittivity  in  the 
field  subject  to  increasingly  high  frequencies,  and  r  is  the  relaxation  time. 


In  case  n  =  2,  we  obtain  the  Lorentz  model 


d2P 
dt 2 


dP 

+  a'~di  +  a2  p  = 


a0E 


or  with  the  physical  parameters 


where  u>q  is  the  natural  frequency  and  uip  is  a  constant  related  to  the  charge  of  the 
medium  oscillator.  In  this  model,  media  are  treated  as  a  fine-grained  assembly  of 
molecular  oscillators.  Such  a  treatment  has  gained  substantial  support  due  to  our 
knowledge  of  molecular  and  atomic  structure  [11]. 

We  study  a  general  n-th  order  ODE  polarization  equation  since  there  is  still  a  sub¬ 
stantial  lack  of  understanding  related  to  the  inability  to  explain  the  electromagnetic 
wave  propagation  phenomena  in  some  media  by  employing  one  of  the  three  dispersion 
models  mentioned  above. 


Next,  we  define  a  one-dimensional  inhomogeneous  dispersive  slab  occupying  [zi,  z2\ 
such  that  0  <  zi  <  z2  <  1  as  our  example  geometry.  The  medium  outside  the  slab 
is  filled  with  air  in  which  e0,  /x0,  <7  =  0.  We  assume  absorbing  boundary  conditions; 


that  is, 


f(t,0)-cf(t,0)  =  0, 
f(M)  +  Cf(U)  =  0, 


where  c  is  the  speed  of  light  in  vacuum. 

It  will  be  shown  that  under  certain  conditions  the  system  (5)  and  (6)  coupled  with 
boundary  conditions  (7)  is  well  posed  in  a  properly  defined  space  in  the  context  of  a 
variational  formulation. 

The  identification  of  electric  polarization  mechanisms  is  formulated  as  parameter  es¬ 
timation  problems.  For  given  incident  (pulse  modulated  microwave  signals,  for  the 
importance  of  such  input  signals,  see  [2])  and  scattered  electric  field,  the  parameter 
estimation  problems  it  to  find  a  set  of  parameters  such  that  the  estimated  electric 
field  matches  the  measured  data  optimally  under  some  criterion.  Let  the  collection  of 
unknown  parameters  be  denoted  by  a  vector  q  =  (a ,  a0, . . . ,  an).  For  given  observa¬ 
tions  {Ei,j}  corresponding  to  measurements  at  times  ft-  and  position  zj,  we  consider 
the  least  squares  estimation  problem  of  minimizing  over  q  €  Q  the  least  squares 
functional  2 

J(E,  E-,q)=  |  {E(U,  Zf,  q)}  -  {Ei,j}  |  ,  (S) 

where  {E(U,  zy,  q)}  are  the  parameter  dependent  weak  solutions  of  (5)-(7)  evaluated 
at  each  time  U,  i  =  1,2,  — ,  and  each  position  Zj,  j  =  1,  •  •  • ,  Nz,  and  |  -  |  is  an 
appropriately  chosen  Euclidian  norm.  The  set  Q  is  some  admissible  parameter  set. 

The  minimization  in  our  parameter  estimation  problems  involves  an  infinite  dimen¬ 
sional  state  space  and  an  infinite  dimensional  admissible  parameter  set  (e.g.,  of  func¬ 
tions  cr(z)).  We  thus  consider  Galerkin  type  approximations  in  the  context  of  the 
variational  formulation  of  (5)-(7).  Solving  the  approximate  estimation  problems,  we 
obtain  a  sequence  of  estimates  {qN,M}.  Parameter  estimate  convergence  and  contin¬ 
uous  dependence  (with  respect  to  the  observations  {Eij})  results  can  be  obtained 
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under  certain  assumptions. 

The  convergence  results  allowing  approximation  of  the  parameter  set  provides  a  sound 
necessary  basis  for  the  reconstruction  of  conductivity  and  coefficients  in  polarization 
equation  (6). 

It  is  well  known  that  evaluation  of  electric  and  magnetic  fields  in  the  microwave 
range  is  computationally  expensive.  To  minimize  the  number  of  function  evalua¬ 
tions,  a  Trust  Region  algorithm  [8]  is  employed  in  solving  the  minimization  problem. 
Moreover,  the  algorithm  has  potential  to  provide  a  global  minimal. 

Our  research  in  electromagnetic  inverse  problems  is  at  its  preliminary  stage.  Our  goal 
is  to  demonstrate  the  ability  to  reconstruct  the  physically  related  parameters  in  the 
time  domain  with  the  expectation  of  (a)  incorporating  the  temperature  dependence  of 
the  structure  into  our  model,  (b)  using  various  type  of  incident  signals  such  as  short 
pulses  and  broad  band  continuous  signal,  and  (c)  utilizing  the  transient  responses  to 
improve  the  quality  of  identification. 
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Abstract 

Object  recognition  can  be  parametrized  systematically  through  physically  robust  wave  objects  by 
linking  features  (observables)  in  scattered  field  data  with  features  on  the  object  (target)  giving 
rise  to  the  data.  The  wave  objects  are  broadly  separated  into  global  (mode)  and  local  (wavefront) 
categories.  Their  parametrization  requires  different  wave-oriented  signal-processing  algorithms 
which  are  implemented  conveniently  in  relevant  subdomains  of  the  configuration  (space-time)  - 
spectrum  (wavenumber-frequency)  phase  space  [1].  We  have  found  that  confidence  in  any  phase- 
space  processing  scheme  is  gained  through  "proper"  physical  understanding  of  scattering 
phenomenologies,  and  that  the  signal-processing  algorithms  can  be  fully  calibrated  by  applying 
them  to  fully  controlled  and  quantified  canonical  prototypes  [2].  We  have  illustrated  these 
processing  strategies  on  several  examples  [1,3,4].  Because  of  space  limitations,  we  discuss  here 
only  the  numerical  techniques  that  have  been  employed  to  generate  the  data  base.  The  examples 
include  short-pulse  scattering  from  a  perfectly  conducting  wire  [3]  or  body  of  revolution  (BOR) 
[4]  buried  in  a  lossy,  dispersive  halfspace  (see  Sec.  I).  The  forward  modeling  is  performed  via 
the  Method  of  Moments  (MoM)  and  allows  for  the  study  of  short-pulse  scattering  from  fairly 
realistic  (though  specially  oriented)  three-dimensional  targets,  for  example,  a  perfectly  conducting 
buried  mine.  In  Sec.  II  we  consider  scattering  from  a  target  in  a  halfspace  with  inhomogeneities 
in  the  soil  profile  and/or  interface  roughness;  to  account  for  such  situations,  we  have  developed 
an  FDTD  code  (so  far  in  two  dimensions).  Concerning  the  inverse  problem,  we  have  performed 
a  systematic  investigation  of  the  resonant  frequencies  of  buried  targets.  The  late-time  resonant 
frequencies  of  targets  in  unbounded  space  are  aspect  independent  and  therefore  have  been  utilized 
for  target  identification.  However,  for  targets  buried  near  the  soil  interface,  the  resonant 
frequencies  will  in  general  be  aspect  dependent  due  to  target-interface  reverberations.  This  issue 
is  investigated  in  detail  in  [3].  We  have  also  performed  space-wavenumber  processing  on  buried- 
target  scattering  data,  with  the  phase-space  signature  calculated  using  the  windowed  Fourier 
transform,  wavelet  transform,  and  windowed  superresolution  processing.  All  of  the  issues  omitted 
here  are  dealt  with  in  the  oral  presentation  or  in  the  cited  references.  Concerning  analytical 
aspects,  the  wave-oriented  modeling  utilizes  local-global  decomposition  (hybrid  ray-mode  and 
wavefront-resonance  schemes),  Gaussian  steady-state  and  pulsed  beams,  deterministic-stochastic 
interactions  [5],  and  others.  The  wave  objects  associated  with  these  models  can  be  forward  and 
backward  propagated  for  task-oriented  reconstruction;  these  aspects  are  now  under  consideration. 
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I.  Short-Pulse  Scattering  from  Buried  Bodies  of  Revolution:  MoM  Solution 


Most  previous  investigations  of  ground  penetrating  radar  (GPR)  phenomenology  have  been 
based  on  narrowband  operating  conditions.  Recently,  however,  there  has  been  an  interest  in  ultra- 
wideband  short-pulse  inputs  for  ground  penetrating  applications,  because  of  the  accompanying 
space-time  resolution.  Such  systems  generate  short  pulse  waveforms  with  large  instantaneous 
bandwidth,  giving  rise  to  time-domain  phenomenology  which  is  fundamentally  different  from  that 
of  narrowband  systems.  As  a  first  step  toward  understanding  the  wave  physics  associated  with 
short-pulse  scattering  from  buried  targets,  we  consider  here  short-pulse  plane-wave  scattering 
from  buried  perfectly  conducting  bodies  of  revolution  (BOR).  To  make  such  a  study  tractable, 
the  soil  is  modeled  as  a  lossy,  dispersive  halfspace,  and  the  axis  of  revolution  of  the  buried  target 
is  assumed  normal  to  the  air-ground  interface;  this  assumption  restricts  the  target  orientation  but 
allows  one  to  view  the  target-halfspace  composite  as  a  single  BOR  with  consequent  simplification 
in  analysis  and  numerics.  We  solve  the  scattering  problem  in  the  frequency  domain  via  MoM, 
using  an  algorithm  analogous  to  those  developed  previously  for  scattering  from  BORs  in 
ffeespace.  After  performing  the  MoM  analysis  over  the  bandwidth  of  the  incident  pulse,  the  time- 
domain  scattered  fields  are  synthesized  via  Fourier  transform. 

Several  authors  have  investigated  the  MoM  analysis  of  plane-wave  scattering  from 
conducting,  dielectric,  and  composite  BORs  in  free  space.  These  algorithms  take  advantage  of 
the  azimuthal  periodicity  inherent  in  the  incident  fields,  induced  currents,  and  the  free-space 
Green’s  function.  Scattering  from  a  buried  BOR  is  complicated  significantly  by  the  fact  that  the 
half-space  Green’s  function  is  a  dyadic,  each  term  of  which  is  expressed  in  terms  of  highly 
oscillatory  Sommerfeld  integrals  and  their  spatial  derivatives.  Thus,  to  make  the  MoM  analysis 
of  a  buried  BOR  tractable,  highly  efficient  techniques  are  required  for  the  computation  of  the 
dyadic  half-space  Green’s  function,  especially  for  the  wideband  applications  here. 

To  compute  the  components  of  the  half-space  Green’s  function  efficiently  over  wide 
bandwidths,  we  employ  the  method  of  complex  images,  developed  by  Chow  et  al.  [3,4].  In  this 
technique,  the  spectral  Green’s  function  inside  the  Sommerfeld  integral  is  fit  via  Prony’s  method 
(or  any  other  parametric  algorithm,  such  as  the  matrix-pencil  method)  to  a  finite  sum  of  complex 
exponentials,  the  integrals  of  which  can  be  expressed  in  closed  form  via  Sommerfeld ’s  identity 
[3,4].  This  formulation  is  highly  efficient  because  the  numerical  burden  is  shifted  from 
laboriously  computing  the  Sommerfeld  integrals  via  brute-force  numerical  integration  to  the 
relatively  efficient  task  of  performing  a  parametric  fit  to  the  spectral  Green’s  function.  For  the 
BOR  problem,  after  applying  the  method  of  complex  images,  each  component  of  the  dyadic 
Green’s  function  is  written  as  a  sum  of  free-space  Green’s  function  (or  its  spatial  derivatives) 
with  complex  source  points.  Thus,  after  applying  the  method  of  complex  images,  the  techniques 
used  for  the  MoM  analysis  of  scattering  from  BORs  situated  in  freespace  -  which  take  advantage 
of  properties  of  the  freespace  Green’s  function  -  can  be  transferred  directly  to  the  problem  of 
scattering  from  buried  BORs. 

For  illustration,  we  consider  short-pulse  scattering  from  a  perfectly  conducting  sphere  and 
cylinder  buried  within  a  lossy,  dispersive  half  space.  The  incident  pulse  is  described  by  the 
waveform  in  Fig.  1,  using  a  free-space  center  wavelength  Xc=l  m  (fc=300  MHz)  consistent  with 
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many  ultra-wideband  radar  systems  of  interest.  The  frequency-dependent  complex  permittivity 
is  taken  from  measured  data  for  Peurto  Rico  clay  with  10%  water  [3,4].  Although  the  algorithm 
is  general,  for  these  examples  we  assume  that  the  pulsed  plane-wave  is  incident  normal  to  the 
air-ground  interface.  Further,  results  are  presented  for  the  scattered  fields  at  the  interface,  since 
it  is  there  that  many  GPR  systems  perform  their  measurements.  We  therefore  do  not  employ  a 
far-zone  approximation;  rather,  we  calculate  the  scattered  fields  using  the  complete  dyadic 
Green’s  function  applied  in  the  MoM  computations  -  again  utilizing  the  complex  image 
technique.  For  airborne  GPR  systems,  one  is  usually  interested  in  the  far-zone  fields,  which  can 
be  calculated  more  easily  by  taking  advantage  of  asymptotic  simplifications  to  the  Green’s 

function. 

The  perfectly  conducting  sphere  has  a  30  cm  diameter,  and  its  center  is  situated  30  cm 
from  the  air-ground  interface.  The  time-domain  fields  for  this  example  are  shown  in  Fig.  2.  For 
the  parameters  considered,  the  first  creeping  wave  response  (also  present  in  the  freespace  case) 
is  clearly  visible,  arriving  after  the  wavefront  scattered  by  the  initial  interaction  of  the  incident 
wave  with  the  target.  For  this  example,  multiple  reverberations  between  the  surface  of  the  sphere 
and  the  air-ground  interface  contribute  negligibly  to  the  total  scattered  field  (we  detect  no 
scattered  fields  at  times  for  which  such  multiple  reverberations  are  expected).  With  the  incident 
wave  from  the  air  described  in  Fig.  1,  the  incident  waveform  that  hits  the  target  is  determined 
by  the  transmission  coefficient  at  the  interface  as  well  as  the  propagation  properties  (dispersion 
and  loss)  between  the  interface  and  the  target. 

The  next  example  considers  short-pulse  scattering  from  a  perfectly  conducting  cylinder 
of  radius  12.8  cm  and  thickness  8.5  cm  buried  inside  the  Peurto  Rican  soil  used  in  Fig.  2;  the 
distance  from  the  top  of  the  cylinder  and  the  air-ground  interface  is  29.8  cm.  The  scattered  fields 
for  this  example,  observed  along  the  axis  of  the  cylinder,  are  shown  in  Fig.  3.  Again,  the 
scattered  fields  are  characterized  by  a  strong  return  from  the  top  of  the  cylinder  followed  closely 
in  time  by  a  weaker  return.  However,  now,  this  later  response  is  attributed  to  a  reverberation 
between  the  air-soil  interface  and  the  flat  surface  at  the  top  of  the  target;  this  additional 
waveform  after  the  specular  return  was  not  visible  in  the  scattered  signature  from  the  same  target 
in  freespace  (with  the  physical  size  of  the  cylinder  modified  such  that  it  was  the  same  size 
electrically  with  respect  to  the  incident  waveform).  For  the  target  orientation  and  excitation 
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Figure  2.  Time-domain  fields  scattered  from  a  perfectly  conducting  sphere  of  radius  15  cm 
buried  in  a  lossy,  dispersive  half  space.  The  center  of  the  sphere  is  30  cm  from  the  soil-air 
interface,  the  angle  of  incidence  is  0-0°,  the  scattered  fields  are  observed  at  the  interface  directly 
above  the  sphere  center,  and  the  polarization  of  the  plotted  scattered  fields  is  the  same  as  that 
of  the  incident  wave.  The  arrow  identifies  the  first  creeping  wave  arrival. 


Figure  3.  Time-domain  fields  scattered  from  a  perfectly  conducting  cylinder  of  12.8  cm  radius 
and  8.5  cm  thickness  buried  in  a  lossy,  dispersive  half  space.  The  top  of  the  cylinder  is  29.8  cm 
from  the  air-soil  interface;  the  angle  of  incidence  is  0;=O°;  the  scattered  fields  are  observed  at  the 
air-soil  interface,  along  the  cylinder  axis;  the  polarization  of  the  plotted  scattered  fields  is  the 
same  as  that  of  the  incident  wave.  The  arrow  identifies  the  first  target-interface  reverberation. 
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II.  Short-Pulse  Scattering  from  Buried  Targets:  FDTD  Solution 

For  this  class  of  problems,  we  have  implemented  a  two-dimensional  FDTD  algorithm 
which  utilizes  a  Huyghens  surface  for  the  incident  field  and  a  second-order  Mur  absorbing 
boundary  condition  (ABC).  The  GPR  problem  involves  a  difficulty  which  is  not  encountered  for 
short-pulse  scattering  from  targets  in  free  space.  Namely,  to  model  plane-wave  incidence,  one 
must  modify  the  FDTD  algorithm  to  simulate  an  infinite  soil  region  without  introducing  artifacts 
at  the  terminations  of  the  model  soil.  This  has  been  done  by  placing  the  model  soil  in  direct 
contact  with  the  ABC.  However,  one  must  now  use  a  Huyghens  surface  for  the  incident  field 
which  simulates  a  pulsed  plane  wave  interacting  with  an  infinite  half  space.  Thus,  the  Huyghens 
currents  account  for  the  wave  which  is  transmitted  into  the  soil  as  well  as  the  reflected  wave 
excited  at  the  air-soil  interface;  if  the  geometry  inside  the  Huyghens  surface  is  a  homogeneous 
half  space,  the  fields  outside  the  Huyghens  surface  vanish.  However,  discontinuities  introduced 
into  the  half  space  (inhomogeneities  in  the  soil  profile,  surface  roughness,  and  targets)  excite 
scattered  fields,  and  these  fields  propagate  outside  the  Huyghens  surface  and  can  therefore  be 
separated  from  the  incident  waveform.  It  is  in  this  manner  that  our  FDTD  GPR  calculations  have 
been  performed  and  an  example  is  shown  here. 

The  problem  geometry  is  shown  in  Fig.  4,  with  a  pulsed  plane  wave  incident  at  45°  (the 
shape  of  the  incident  pulse  is  as  in  Figs.  2  and  3).  This  geometry  has  been  selected  to  simulate 
a  rough  soil  surface  and  a  buried  perfectly  conducting  target  located  in  a  region  of  recently 
replaced  excavated  soil  (rectangular  enclosure);  it  has  been  shown  in  numerous  GPR 
measurements  that  recently  replaced  excavated  soil  has  an  index  of  refraction  different  from  that 
of  undisturbed  soil.  The  results  in  Figs.  5  (a)-(c)  demonstrate  the  space-time  evolution  of  the 
scattered  fields.  We  can  see  that  the  surface  roughness  introduces  interesting  focusing  effects 
inside  the  soil.  Further,  the  scattered  field  from  the  target  in  Fig.  5(c)  is  seen  to  be  affected  by 
phenomena  due  to  the  contrast  in  dielectric  constant.  Such  phenomenological  studies  yield  insight 
into  the  GPR  scattering  and  propagation  mechanisms,  and  therefore  are  useful  in  guiding  the 
development  of  wave-parametrized  signal  processing  and  imaging  architectures. 
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Figure  4.  Geometry  used  in  FDTD  example. 
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Figure  5(a) 


(b)  (c) 

Figure  5.  Space  time  evolution  of  fields  scattered  from  rough-soil-surface  and  buried  target.  Time 
increases  from  (a)  to  (c),  with  each  plot  representing  the  spatial  distribution  of  the  fields  at  a 
snapshot  in  time. 
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1  Introduction. 


In  this  study,  the  maximum  likelihood  estimation  procedure  is  used  to  locate  buried 
pipes  based  on  experimental  ground  penetrating  radar  (GPR)  return  data.  The  report 
begins  with  a  review  of  the  relevent  estimation  theory,  then  describes  how  the  theory  is 
applied  to  our  particular  application.  Finally,  results  are  reported  from  an  experimental 
GPR  survey. 


2  Estimating  parameters  from  noisy  measurements 
-  probabilty  densities  and  log-likelihood  functions. 


This  report  is  based  on  the  general  problem  of  estimating  a  set  of  parameters  from  noisy 
observations.  Here,  the  parameters  are  the  location  coordinates  of  a  buried  object,  while 
the  observations  are  voltage  samples  of  radar  return  pulses.  Assume  the  following  model 
for  the  measurement  of  data  s  in  the  presense  of  additive  noise  n: 

s  =  d(x)  +  n.  (1) 
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Here,  x  =  (xi,  x 2 •,  x  3 ,  ••••  x  n )  is  the  vector  of  parameters  to  estimate  given  the  measured 
data  vectors  =  ($1,  S2,  S3, sm),  while  d(x)  =  (di(x),d2(x),...,dm(x))  is  the  parameter- 
dependent  data  model,  i.e.  the  data  that  would  be  expected  in  the  absence  of  noise. 


A  natural  way  to  estimate  x  would  be  to  construct  the  conditional  probability  density 
function  p(x/s),  to  evaluate  it  for  the  measured  data  s,  then  to  maximize  the  function 
over  the  domain  of  x.  The  function  is  relatively  easy  to  compute  assuming  uncorrelated, 
Gaussian  white  noise.  Then  [1] 


p(x/s) 


P(s/x)p(x) 

P( s) 


(2) 


where,  for  independent  data  elements  s,-, 

m  j 

P( s/x)  =  II 

i= 1  \/2Trcrf 

Here,  the  noise  component  rii  has  mean  p,-  and  variance  <7?. 


-1 
2a ? 


| (si  -  m)  -  di(x)\J 


(3) 


Often,  however,  the  a  priori  probabilities  p(x)  and  p(s)  will  not  be  known.  In 
55 maximum  likelihood  estimation” ,  these  probabilities  are  set  to  unity,  and  x  is  estimated 
simply  by  maximizing  p( s/x)  as  given  in  equation(3).  Or,  since  the  natural  log  function 
is  monotonic,  the  estimate  can  be  made  by  maximizing  the  log  of  equation(3),  i.e. 


m  _  I 

L(*)  =  £  1T2  l(S«  “  Pi)  ~  ^(*)l  • 

j=l  Z<Ji 

This  function  is  known  as  the  log-likelihood  [1]. 


(4) 


An  alternate  expression  results  from  expanding  the  squared  term  and  retaining  only 
the  terms  which  depend  on  x.  Then 

m  1 

m=T,  4(2Si<s<  -  «)•*(*)]  -  &(*)!).  (5> 

t=l  ai 

where  the  *  indicates  the  complex  conjugate  of  any  non-real  values. 


The  log-likelihood  function  L(x)  can  be  computed  over  a  range  of  parameters  x,  i.e. 
a  ”  log-likelihood  image”  can  be  constructed,  where  the  estimate  x  will  correspond  to  the 
location  of  the  image  peak.  If  the  signal-to-noise  ratio  of  the  data  is  high  enough,  the 
peak  of  the  log-likelihood  image  should  roughly  correspond  to  the  actual  object  location. 
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Note  that  the  first  term  in  equation(5)  is  simply  the  correlation  between  the  measured 
data  and  the  data  model  (after  subtracting  /*,•  and  weighting  by  jj).  The  second  term, 
which  may  be  called  the  ’’energy  term” ,  is  often  slowly  varying  over  the  domain  x  relative 
to  the  correlation,  and  therefore  may  not  significantly  affect  the  maximum  likelihood 
peak  within  a  local  radius.  Therefore,  in  certain  cases  the  log-likelihood  function  is 
simply  the  weighted  correlation  between  the  measured  data  and  the  data  model. 


3  Estimation  theory  for  GPR  surveys. 


The  general  notation  of  the  preceding  section  is  now  adapted  to  the  particular  application 
of  time  domain  GPR.  Suppose  that  a  GPR  survey  consists  of  a  series  of  experiments, 
from  which  it  is  hoped  to  locate  a  buried  target  object  with  known  characteristics  (e.g. 
a  six  inch  steel  pipe).  An  "experiment”  is  defined  as  a  collection  of  measurements  made 
at  a  single  receiver  location  simultaneous  with  the  excitation  of  a  single  transmitter. 
The  experiments  can  be  either  monostatic  or  bistatic.  At  each  receiver,  measured  data 
consists  of  a  time  domain  pulse,  scattered  by  the  target  object,  corrupted  by  noise,  and 
then  sampled.  Thus,  the  following  signal  model  results  in  the  time  domain: 

sjk(ti)  =  djk(ti;  x0)  +  njk(ti),  (6) 

where 

•  Sjk(ti)  is  the  noisy  signal,  sampled  at  time  th  and  measured  at  receiver  k  simulta¬ 
neous  to  the  excitation  of  transmitter  j; 

•  djk(ti',x o)  is  the  data  model  for  an  object  located  at  coordinates  Xo  =  (^O;  Vo-t  zo)] 

•  rijk(ti)  is  additive  Gaussian  white  noise  with  variance  <r]k(ti)  and  mean 

As  described  above,  the  log-likelihood  function  has  the  following  expression: 

L(x)  =  £  £  £  -sVt  [2M<-)  -  *)  -  <?*(«■;  *)]  •  (7) 

j= 1  k= i  1=1  VjkW) 

It  is  evident  that  the  total  log-likelihood  function  can  be  computed  by  adding  the  indi¬ 
vidual  log-likelihood  functions  from  each  of  the  J xK  experiments. 
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4  A  simple  data  model  for  time  domain,  monostatic 
experiments. 


In  the  time  domain,  a  simple  data  model  can  be  derived  which  is  valid  for  monostatic 
experiments,  and  is  based  on  the  following  approximations: 

•  the  soil  and  buried  object  are  nondispersive  (frequency  independent); 

•  the  target  object  is  far  enough  from  the  antennas  that  radiation  incident  on  the 
object  looks  locally  like  a  plane  wave,  as  does  radiation  scattered  to  the  receiving 
antenna; 

•  the  target  object  is  small  compared  with  all  wavelengths  of  the  radar,  thus  it  can 
be  treated  as  a  point  scatterer; 

•  the  radiation  pattern  of  the  antennas  is  independent  of  the  radar  frequency. 

Using  the  approximations,  for  a  receiver  located  at  the  origin  the  data  model  can  be 
expressed  as 

d(t- x)  =  A\v )p(t  -  2 r/c)F2(r)D,  (8) 


where 


•  r  is  the  length  of  the  vector  x,  while  f  is  a  unit  vector  with  the  direction  of  x; 

•  2 r/c  is  the  round-trip  propagation  time  for  the  pulse  (c  is  the  propagation  velocity); 

•  p(t)  is  the  second  derivative  of  the  transmitted  pulse  [5]; 

•  F(r)  is  a  radial  gain  function  consisting  of  geometrical  spreading,  soil  attenuation, 
and  any  other  losses  (e.g.  energy  loss  from  ground  surface  reflection); 

•  D  is  a  quantity  representing  the  scattering  strength  of  the  target; 

•  A(r)  gives  the  directional  dependence  of  the  power  radiated  by  the  antennas. 
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Using  equation(7),  the  time  domain  log-likelihood  function  for  each  experiment  is 
then 


L(x)  =  A2(t)F2(t)D  £  -  K‘i)W i  -  2r/c),  (9) 

1=1  a  uu 


where  the  energy  term  has  been  ignored.  Notice  that  the  directional  r  dependence  of  this 
expression  is  a  multiplicative  factor.  Thus  the  log-likelihood  function  for  a  single  exper¬ 
iment  can  be  quickly  computed  by  first  evaluating  L  along  a  radius  in  a  single  direction, 
then  sweeping  this  one  dimensional  vector  over  all  directions  f,  and  finally  multiplying 
by  the  directionally  dependent  function  A(r).  The  total  log-likelihood  function  is  then 
obtained  by  adding  the  subimages  from  all  monostatic  experiments. 


4.1  Estimating  the  noise  variance  and  radial  gain  function  from 
measured  data. 

In  preceding  section  a  method  was  outlined  for  computing  log-likelihood  functions  in 
time  domain  GPR.  However,  the  method  depends  on  knowing  both  the  noise  statistics, 
c r2(tk )  and  and  the  soil-dependent  gain  function  F(r).  Due  to  the  variable  nature 

of  soil  (electrical  properties  and  clutter  content),  these  functions  cannot  be  accurately 
predicted  in  practice.  However,  we  can  estimate  them  from  our  data. 

Based  on  observations  of  actual  GPR  return  signals,  it  is  evident  that  noise  results 
mainly  from  ground  clutter.  For  a  particular  experiment  in  a  GPR  survey,  the  return 
data  will  be  mostly  noise  if  the  target  object  is  in  the  periphery  of  the  antenna.  This 
will  be  the  case  with  most  experiments  in  a  typical  survey,  therefore  the  noise  mean  and 
variance  can  be  estimated  simply  by  averaging  across  all  experiments  in  a  survey. 

It  is  also  possible  to  approximate  the  radial  gain  function  F(r )  by  using  the  estimated 
noise  variance.  Consider  the  subsurface  of  the  earth  as  a  distribution  of  tiny  clutter 
volumes,  where  scattering  from  each  volume  contributes  to  the  total  noise.  Assume 
that  the  scattering  potential  of  each  volume  is  a  zero-mean  random  process,  which  is 
reasonable  if  air  pockets  as  well  as  rocks  may  be  embedded  in  the  soil.  The  variance  of 
the  noise  component  from  each  tiny  volume  is  governed  by  a  relation  similar  to  that  of 
the  target  signal  in  equation(8),  i.e., 

d„(t-,X)  =  \A2(t)F\ctft)]2*l  (10) 
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where  the  transmitted  signal  is  assumed  to  be  a  very  short  pulse  p(t)  —  6 (t).  Here,  <rc 
is  the  variance  of  the  scattering  potential  for  each  tiny  volume  of  clutter. 


For  a  short  pulse,  at  each  time  t,  the  total  noise  received  originates  from  clutter 
contained  in  a  cylindrical  shell  with  radius  r  =  cf/2,  i.e.  with  volume  proportional  to 
r  =  (d/2)  (we  are  treating  the  problem  as  2D).  Therefore,  the  total  noise  variance  at 

time  t  is  described  by  the  relation 

ol(t)  cc  (ct/2)[F2(ct/2)]2.  (11) 


The  radial  dependence  of  the  gain  function  is  therefore  given  by  the  relation 


F2(r)  oc 


o-n 


(2r/c) 

y/r 


(12) 


It  should  be  noted  that,  in  practice,  this  function  will  diverge  as  r  goes  to  zero.  Since 
it  was  derived  using  some  farfield  approximations,  it  should  therefore  not  be  used  for 
small  radial  distances.  Work  is  being  done  on  near-field  modeling. 


5  Experimental  results. 

A  multimonostatic  GPR  survey  was  conducted  using  a  ’’Georadar”  step  frequency  radar 
with  a  frequency  range  of  450  steps  from  100MHz  to  1GHz  (2MHz  increments).  The 
test  site  consisted  of  fairly  wet  soil  in  which  metal  pipes  were  buried  at  depths  of  less 
than  two  meters.  The  survey  was  performed  by  moving  the  radar  set  continuously  in 
straight  lines  perpendicular  to  the  axis  of  the  pipes  (the  direction  of  motion  corresponds 
to  the  x  axis  on  all  image  plots).  Each  experiment  (record)  in  the  survey  corresponds 
to  data  measured  at  a  certain  x  coordinate.  Radar  return  data  thus  consists  of  450 
complex  amplitudes  (magnitude  and  phase  at  each  frequency)  per  experiment.  In  order 
to  process  in  the  time  domain,  the  complex  amplitude  information  (frequency  domain) 
was  preprocessed  by  Fourier  transforming. 

Results  from  the  various  stages  of  processing  are  shown  in  Fig.l  for  a  particular 
sub-survey.  There  are  known  to  be  two  pipes  in  the  corresponding  section  of  ground,  at 
coordinates  (*,*)  =  (2.8, -0.6)  (4”  diameter)  and  (1.8, -1.0)  (6”  diameter).  The  pipes 
are  evident  in  all  images  of  Fig.  L  The  survey  consists  of  90  experiments  (x  positions), 
with  70  time  domain  return  samples  per  experiment. 
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In  la.,  the  return  pulses  are  simply  stacked  without  range  gaining.  A  surface  reflec¬ 
tion  is  evident  at  shallow  depths  across  all  experiments,  while  the  pipe  at  (2.8,  —0.6) 
appears  as  a  typical  hyperbolic  arc.  It  is  difficult  to  observe  the  second  pipe.  In  lb., 
the  returns  have  been  weighted  by  the  noise  standard  deviation,  averaged  across  all 
experiments  for  each  time  sample.  This  procedure  is  consistent  with  the  discussion  in 
the  preceding  sections.  Now  it  is  much  easier  to  observe  the  second  pipe,  which  ap¬ 
pears  faintly  hyperbolic,  and  seems  to  be  resonating.  In  lc.,  the  log-likelihood  function 
has  been  computed  as  described  by  equations(9)  and  (12),  where  the  transmitted  pulse 
has  been  approximated  as  a  delta  function.  The  directional  dependence  of  the  anten¬ 
nas  has  been  modeled  as  cos2(0),  which  seems  consistent  with  the  data.  The  effect  of 
this  processing  seems  to  be  a  focussing  of  the  hyperbolic  arcs.  Notice  that  the  surface 
reflections  appear  greatly  diminished.  Finally,  in  Id.  the  log-likelihood  function  has 
been  thresholded  at  a  level  which  is  80  percent  of  the  maximum  absolute  value.  This 
result  demonstrates  the  potential  value  of  the  algorithm  for  an  automatic  data  process¬ 
ing  scheme.  Both  pipes  appear  in  the  image,  with  few  false  readings.  Recall  that  the 
processing  has  required  no  operator  adjustments  or  interventions. 
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Abstract 

The  frequency  dependent  electrokinetic  properties  of  rocks  and  soils  are  currently 
unknown,  even  though  they  are  vital  to  our  understanding  of  the  electric  signals  that 
are  observed  to  be  generated  by  seismic  waves  propagating  in  the  earth.  These  electric 
fields  are  caused  by  relative  fluid  flow  in  rocks,  induced  by  a  passing  seismic  wave. 
This  relative  flow  then  induces  a  streaming  current  that  results  in  an  electric  dipole, 
which  is  the  source  of  a  electromagnetic  wave  that  oscillates  with  the  same  frequency 
as  the  seismic  excitation. 


Electroseismic  signals  have  been  measured  in  a  variety  of  field  environments,  and 
corresponding  theoretical  studies  are  currently  under  way.  The  differences  in  the 
electroseismic  prospecting  (ESP)  and  the  seismic  sections  can  reflect  the  fact  that  the 
ESP  signals  are  not  only  sensitive  to  the  seismic  structure,  but  are  also  sensitive  to  the 
electrokinetic  properties  of  the  formation  (Thompson  and  Gist,  1991,  1993). 
Electroseismic  prospecting  is  potentially  a  very  powerful  new  exploration  method 
that  may  provide  information  both  about  subsurface  hydrogeological  parameters  (e.g., 
porosity  and  permeability)  and  about  the  subsurface  chemical  environment. 

The  coupling  of  hydraulic  flow  and  electric  currents  within  the  earth  is  governed  by 
Onsager’s  laws  (1931)  of  irreversible  thermodynamics,  namely: 


q  =  —  VP  —  c  Cv  VV 
r\ 

(1) 

)  =  -arCvVP  -  cr  VV , 

(2) 

where  q  is  Darcy’s  volume  fluid  velocity  per  unit  cross-sectional  area  and  j  is  the 

electrical  current  density.  VP  and  VV  are  the  hydraulic  pressure  gradient  and  the 
electrical  potential  gradient,  receptively.  The  rock  permeability  is  given  by  k  and  its 

electrical  conductivity  by  ar.  The  viscosity  of  the  pore  liquid  is  t|.  Cv  is  the  voltage 
cross-coupling  coefficient  and  gives  the  voltage  generated  per  unit  pressure  gradient. 
Fluid  flow  problems  within  the  earth  can  be  modeled  by  self-consistently  solving 
equations  (1)  and  (2)  with  the  appropriate  boundary  conditions  (Nourbehecht,  1963; 
Sill,  1983;  Wurmstich  etal.,  1991;  Wurmstich  and  Morgan,  1992,  1994).  The  major 
weak  point  in  solving  equations  (1)  and  (2),  even  for  steady  state  flows,  is  the  sparcity 
of  Cv  data. 

To  better  understand  and  model  the  electroseismic  phenomena,  we  must  know  the 
frequency  dependence  of  the  electrokinetic  coupling  coefficients  for  various  rock  and 
soil  types.  The  limited  data  that  exist  in  the  chemistry  literature  for  capillaries  and 
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porous  diaphragms  indicate  that  the  streaming  potential  is  very  sensitive  to  frequency 
and  that  it  is  a  complex  function  of  the  capillary  geometries  and  the  solution 

chemistry.  . 

Packard  (1953)  measured  the  streaming  potential  responses  of  soda  glass  capillaries 
of  varying  radii  (1.67  mm  and  1.46  mm)  (figure  la).  Packard  developed  a  theory  for 
the  frequency  dependent  streaming  potential  responses  of  capillaries.  His  theory 

predicts  that  for  frequencies  greater  than  a  critical  pulsation  coo  the  streaming 
potential  response  will  decrease  with  increasing  frequency  (figure  lb).  For 
frequencies  less  than  too,  the  streaming  potential  response  asymptotically  approaches 
the  (dc)  Helmholtz-Smoluchowski  value.  The  critical  frequency  is  given  by 


where  t)  and  p  are  the  fluid  viscosity  and  density,  and  a  is  the  capillary  radius.  The 
critical  frequencies  for  the  G2  and  G4  capillaries  are  0.051  Hz  and  0.067  Hz, 
respectively. 

Cooke  (1955)  measured  the  streaming  potential  responses  of  porous  diaphragms  over 
the  frequency  range  of  1  Hz  to  400  Hz  (figure  2).  The  samples  were  saturated  with 
0.01  N  KC1  solutions.  The  critical  frequency  is  approximately  20  Hz,  which  is  much 
larger  than  the  critical  frequencies  for  the  capillaries. 

Streaming  potential  measurements  were  made  by  Sears  and  Groves  (1977)  on 
capillaries  coated  with  bovine  serum  albumin  (BSA)  (figure  3).  Measurements  on  a 
capillary  with  a  radius  of  a  =  0.508  mm  showed  a  critical  frequency  of  0.56  Hz.  The 
streaming  potential  data  begin  to  flatten  off  at  a  slightly  higher  frequency  of 
approximately  1 .0  Hz. 

This  frequency  dependence  should  be  even  more  pronounced  in  porous  media,  such 
as  rocks,  in  which  there  are  wide  distributions  of  pore/grain  sizes.  Recently  the 
governing  equations  for  the  coupled  electromagnetics  and  acoustics  of  porous  media 
has  been  derived  (Pride,  1995)  and  the  critical  pulsation  of  the  streaming  potential  in 
a  porous  medium  has  been  shown  to  be 

©o  =  ti  /  (F  ko  pf )  (4) 

with  F  the  formation  factor  measured  at  high  frequency,  k0  the  permeability 

measured  at  zero  frequency  and  pf  the  density  of  the  fluid.  This  critical  frequency  is 
the  same  as  the  critical  frequency  in  the  dynamic  permeability  theory  (Johnson  et  al., 
1987).  Knowing  that  we  can  expect  higher  critical  frequencies  as  suggested  by  some 
dynamic  permeabilty  measurements  on  glass  beads,  sandgrains  (Smeulders  et  al., 
1992)  or  on  fused  glass  beads  and  crushed  glass  (Charlaix  et  ah,  1988).  These  authors 
measured  critical  frequencies  of  the  order  of  100  Hz  for  very  high  permeable 
synthetic  samples,  and  we  expect  higher  critical  frequencies  for  rocks  and  soils. 

Briefly,  we  propose  a  project  with  the  following  three  parts:  (1)  to  design  and 
construct  an  experimental  system  to  measure  frequency-dependent  electrokinetic 
properties  of  rocks  and  soils;  (2)  to  measure  the  cross-coupling  coefficients  of  rock 
and  soil  samples  as  a  function  of  frequency  from  dc  through  the  seismic  frequency 
band  (1  kHz),  sample  microgeometry,  and  solution  chemistry;  and  (3)  to  develop  a 
theory  for  the  frequency  dependence  of  the  electrokinetic  properties  of  rocks  and 
soils.  This  research  will  be  crucial  in  improving  our  understanding,  modeling,  and 
interpretation  of  the  electroseismic  phenomena  and  in  determining  its  utility  in 
applications,  such  as  hydrological  characterization,  contaminant  detection,  and 
formation  evaluation  within  a  borehole  environment. 
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Figure  la:  Effective  value  of  streaming  potential  versus  frequency  for  a  constant 
pressure  amplitude  of  10,000  Dynes/cm2  (0.1  atm)  (from  Packard,  1953). 


Figure  lb:  Comparison  of  experimental  values  of  streaming  potential  versus 
frequency  with  theory.  Ya  is  a  reduced  frequency  variable,  where  co  is  the  angular 
frequency,  a  is  the  capillary  radius,  p  is  the  fluid  density,  and  p  is  the  fluid  viscosity 
(from  Packard,  1953). 


Figure  2:  Plots  of  streaming  potential  versus  frequency  for  three  porous  diaphragms 
(from  Cooke,  1955) 


Figure  3:  Streaming  potential  (cross-coupling  coefficient)  versus  frequency  for  a  BSA- 
coated  capillary  (from  Sears  and  Groves,  1977). 
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Abstract 

In  this  abstract,  we  are  concerned  with  characterizing  the  internal  structure  of  2D  strongly 
scattering  objects,  using  a  single  illumination  direction  or  view  and  constant  wavenumber 
k0=2%  /  L0 .  The  object  is  illuminated  with  the  plane  wave  e'k°TrT  and  simulated  scattered 
far-f  ds  are  calculated  on  a  circular  aperture  in  the  plane  of  the  object's  cross  section.  A 
differential  cepstral  filter  is  applied  to  a  single-view  backpropagated  image  which  relates 
to  the  product  of  the  scattering  potential  and  the  total  field.  This  nonlinear  filter  avoids 
the  phase  wrapping  problems  associated  with  homomorphic  filtering,  and  is  used  to 
isolate  the  scattering  potential  from  single-view  backpropagated  images. 

The  scalar  Green's  function  solution  to  the  scattered  field  from  an  inhomogeneous  region 
can  be  written  1 


¥'(!■„?,)  =  -k]  jF(r')'F(r',r,)G(rf  ,r')*'  (1) 

D 

where  V(r)  is  the  scattering  potential,  ¥(r,r,.)  is  the  total  (incident  ¥'(!%?,)  plus 
scattered  '¥s(r,rl))  field,  G(rs,r)  is  the  free  space  Green's  function,  r,  indicates  the 
illumination  direction,  and  the  integral  is  over  the  object's  compact  support  D.  In  two 

dimensions,  G(rs =  |rf  - r|)  where  J^)(A:0|rJ  —  r|)  represents  a  zero  order 

Hankel  function  of  the  first  kind.  When  |r5  -  r| »  Z)  (far-field)  one  may  use  the 
following  approximation 


^o!^oK-r|) 


4 

i 


I  1 

V  k0rs 


(2) 


In  this  case,  the  scattered  far-field  in  equation  (1)  may  be  written 
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n<r,,r,.)  =  - 


(3) 


l2kKi(^rs,r,) 


where  the  scattering  amplitude  F(rJ,r/)  is  defined  as 


(4) 


It  is  common  to  adopt  weakly  scattering  approximations  to  solve  equation  (4).  For 
example,  under  the  Bom  approximation  (BA),  one  uses  the  incident  field  as  an 
approximation  to  the  total  field  reducing  the  scattering  amplitude  to  the  following^ 


Fba(  rs,r,)  =  — — y  J  V(r')e-iioCr-f'>r'dr' . 
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(5) 


The  Fourier  transform  of  the  scattering  potential  is  defined  as 


F(k)  =  — J-y  [  V(r')e~ikr  dr' 
4  n  J 


(6) 


and  we  recognize  that  FfM(rJ,r/)  specifies  the  Fourier  transform  of  the  scattering 
potential  at  k  =  k0(rs  -r,).  An  image  of  V(r)  may  be  reconstructed  by  either  Fourier- 
based  interpolation  or  by  a  backpropagation  process3. 

In  practice,  most  objects  are  not  weakly  scattering  and  the  Bom  approximation  is  not  a 
valid  one.  For  these  more  strongly  scattering  objects,  the  total  field  T'(r,r,)  in  equation 
(4)  is  difficult  to  determine  a  priori  due  to  internal  multiple  scattering  effects.  This 
makes  the  direct  problem  of  calculating  the  scattered  field  a  difficult  one  and  greatly 
complicates  the  inverse  problem.  However,  by  multiplying  the  integrand  in  equation  (4) 
by  unity  (i.e.  'F'(r,r,.)/ vF'  (r,r,  ))  one  may  formulate  the  scattering  amplitude  as^ 


where. 


F(  = 


1 

47I2 


F(r)Y(r,ff) 

T"(r,f,.) 


(7) 


In  equation  (7),  the  scattering  amplitude  determines  Fourier  data  on  for  those 

frequencies  ^0(fs-f;).  Under  the  Bom  approximation,  F^(r,f,)  reduces  to  F(r)  and 
Fourier  inversion  of  the  scattering  amplitude  data  collected  over  many  illumination 
directions  yields  a  low  pass  estimate  of  V(r).  For  strongly  scattering  objects  we 
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recognize  that  Fourier  inversion  of  equation  (7)  is  not  strictly  valid^  since  Fe#(r,f,  ) 
varies  with  r(.  However,  by  choosing  r,  to  be  constant,  Fourier  data  limited  to  the 
circumference  of  an  Ewald  circle  may  be  acquired. 

The  reconstructed  image,  using  a  single-view  in  the  backpropagation  algorithm^,  may  be 
written 


S(  r'.r,  =  constant)  =  ^J  F(r5,r,  y  W.-W**  d$s 


(8) 


where 


S(r,r ■  =  constant)  = 


r(r)y(r,ff) 

nr,r,) 


!*/?(r,f,) 


and  **  denotes  2D  convolution  with  the  point  spread  function  (PSF)  /2(r,r,).  In 
frequency  space,  the  Fourier  transform  of  h( r,r,.)  can  be  considered  a  window  function 

with  a  value  ^/l  -  (?f  •  r;.  )2  on  the  Ewald  circle  k0(rs  -  r()  and  zero  otherwise.  We  consider 
objects  that  have  inhomogeneities  large  compare  to  the  spatial  resolution  of  /?(r,f,)  so 
that  the  effects  of  the  PSF  can  be  neglected. 


To  isolate  V(r)  in  equation  (8)  from  the  perturbing  factor  vF(r,r,  )/vF'(r,r,  ),  we  utilize  a 
filter  in  the  spectrum  of  the  derivative  of  the  logarithm  of  Sy,?,)  (i.e.  the  differential 
cepstrum).  Direct  linear  filtering  is  not  appropriate  for  multiplied  signals  of  this  kind 
since  their  spectra  are  convolved.  Phase  wrapping  problems  are  avoided  in  the 
differential  cepstrum  by  omitting  the  logarithm  explicitly  in  its  definition.  Neglecting 
h(r,rt),  the  differential  cepstrum7  of  equation  (8)  becomes  the  spectrum  of 


OX _ 

S(  r,r() 


d_V(  .  a  T(r,f;) 

ar  W 

V(r)  Whs) 

'v'M) 


(9) 


After  filtering,  the  reconstructed  image  of  V(r)  is  found  by  integration,  and 
exponentiation.  Successful  recovery  of  F(r)  depends  on  how  well  the  two  terms  on  the 
right  side  of  equation  (9)  separate  in  the  differential  cepstral  domain. 

Consider  two  lossless  concentric  cylinders  with  inner  er  =2,  radius  3.6cm  and  outer 
sr  =  4,  radius  9.9cm,  illuminated  by  a  plane  wave  at  10GHz.  Exact  scattered  far-field 
data  were  calculated  and  backpropagated  (using  equation  (8))  to  a  64x64cm7  physical 
space.  The  magnitude  of  the  single-view  backpropagated  image  is  shown  in  Figure  1 . 
The  scattering  potential  V(r)  is  difficult  to  recognize  since  the  total  field  ^(r,?,.)  is  very 
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different  from  the  incident  field  vP'(r,r/)  in  this  case  (i.e.  object  represents  a  strong 
scatterer). 

Using  a  priori  knowledge  of  the  object's  size,  the  quantity  S(r,r;.),  of  Figure  1,  was  set  to 
zero  just  beyond  its  physical  support  D.  A  simple  low  pass  filter  was  then  applied  to  a 
regularized  differential  cepstrum.  Figure  2  displays  the  reconstruction  of  the  scattering 
potential  generated  from  a  filtered  (in  the  differential  cepstral  domain)  single-view 
backpropagated  image  for  the  same  field-of-view  as  in  Figure  1.  Although  uniqueness  is 
not  necessarily  guaranteed^,  the  reconstruction  qualitatively  agrees  with  the  original 
scattering  potential  V(r). 

In  summary,  single-view  backpropagated  images  are  formulated  as  the  product  of  the 
scattering  potential  and  the  total  field.  A  nonlinear  processing  technique  is  applied  to 
remove  the  perturbing  factor  that  tends  to  complicate  images  of  more  strongly  scattering 
objects.  Limitations  of  the  proposed  technique  arise  as  the  effects  of  the  PSF  become 
more  prominent.  In  such  cases,  the  backpropagated  image  of  r+M  =  constant)  is 
blurred  to  some  extent  due  to  the  convolution  with  h(r,rt  =  constant).  Currently,  we  are 
investigating  a  frequency  extrapolation  routine  that  exploits  a  priori  knowledge  of  the 
object's  support  in  an  effort  to  improve  the  backpropagated  images  from  strongly 
scattering  structures. 


Figure  1.  Figure  2. 
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